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ABSTRACT OF THE DISSERTATION 

Nonlinear Damping Model 
for Flexible Structures 

by 

Weijian Zhang 

Doctor of Philosophy in Electrical Engineering 
University of California, Los Angeles, 1990 

Professor A. V. Balakrishnan, Chair 

This dissertation is on the study of nonlinear damping problem of flexible struc- 
tures. Both passive and active damping, both finite dimensional and infinite dimen- 
sional models are studied. 

In the first part of this dissertation, the spectral density and the correlation 
function of the following single DOF nonlinear damping model is investigated 

x + 2 £u> 0 i + 7 D(x, i) u>qX = an(t) 

where 7 > 0 is a small parameter. A formula for the spectral density is established 
with 0{ 7 2 ) accuracy based upon Fokker-Planck technique and perturbation. The 
spectral density depends upon certain first order statistics which could be obtained 
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if the stationary density is known. A method is proposed to find the approximate 
stationary density explicitly. 

In the second part of this dissertation, the spectral density of the following multi- 
DOF nonlinear damping model is investigated 

Mx + DqX + 7 D(x, x) + Kx = crn(t) 
where 7 > 0 is a small parameter. 

In multi-DOF case, x and x are generally not uncorrelated in stationary state, 
even in linear case, which is one of the features of the multi-DOF model. A necessary 
and sufficient condition for uncorrelatedness is given for the linear model. 

In the third part of this dissertation, energy type nonliner damping model in an 
infinite dimensional setting is studied. According to its geometry of the structures 
considered, the nonlinear damping models are divided into two types. The existence 
and uniqueness result of the nonlinear damping model is based on the work of A. 
Lunardi. 

Then a Krylov- Bogoliubov type approximation is established for the nonlinear 
damping model in the case the linear damping part is neglected. In general, the 
generalization of Krylov- Bogoliubov approximation method, which applies only to 
single DOF model, to multi-DOF model has been a formidable task. The result 
presented here is based upon the specific form of nonlinearity - energy type damping. 
From its Krylov-Bogoliubov approximation, we can see that there is no exchange of 
energy between modes, i.e., internal resonance does not exist. 

The notions of Characteristic equation and its Root locus are extended to actively 
damped distributed parameter systems. The root locus provides an insight of the 
nature of active damping. Sufficient conditions of strong stabilizability are provided, 
which are the weakest sufficient conditions obtained so far. 
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Chapter 1 
Introduction 

1.1 Background and Survey 

A new generation of spacecrafts often contain very large flexible components such 
as truss structures, solar panels, dish antennas, radar arrays, space telescopes and 
space stations. In general, these structures are characterized by weak damping, and 
interconnection of rigid and flexible parts. The tasks of controlling the rotation, 
pointing with high accuracy in minimum time or with minimum energy consump- 
tion, and at the same time stabilizing the vibrations, pose very difficult control 
problems. This need is evident for various ongoing government programs such as 
space shuttle and space station. In 1984, NASA Langley Research Center and Dr. 
A. V. Balakrishnan initiated a design challenge for the SCOLE (Spacecraft Control 
Laboratory Experiment) problem[4], the objective of which includes the task of di- 
recting the line-of-sight of the shuttle/antenna configuration toward a fixed target, 
under condition of noisy data, limited control authority and random disturbances 

[4] - [13]- 
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In most applications, damping is important to the structured dynamics, and in 
many applications it is in fact critical. Two examples may suffice to illustrate the 
point: for spinning (or partly spinning) spacecraft, the level of energy dissipation 
in the structure determine whether an initial wobbling motion grows or decays 
(dynamical stability); and an automatically controlled flexible spacecraft may act 
either unmanageable or docile (control system stability) depending on the level of 
damping in the higher order vibration modes. 

In addition, the advance of modern material science and technology has provided 
us with useful structural material which are generally light weight. Their applica- 
tions in spacrcrafts, high-performance helicopters have sharply increased in the past 
decade. However, such viscoelastic materials have highly nonlinear characteristics 
that cause significant nonlinear response in the system. The questions of analysis, 
design and control appear more difficult. 

Literature Survey 

The following literature survey is made along the lines of nonlinear damping 
and linear damping. In the nonlinear damping case, literature is devided into two 
groups: (1) Finite dimensional model; (2) Infinite dimensional model. While in the 
linear damping case, literature on distributed parameter system is divided into three 
groups: (1) Strictly proportional and asymptotic proportional damping operator; 
(2)Boundary damping model; (3)Interior point damping model. In the end, the 
work on finite dimensional linear damping modelling is briefly reported . 
Nonlinear Damping Model 

• Finite Dimensional Model 

There has been large amount of literature on finite dimensional nonlinear 
damping model, among the authors, T. K. Caughey [25] - [27], S. H. Crandall 
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[34] - [36], P -T. Spanos [64] [65], V. V. Bolotin [21], T. S. Atalik and S. 
Utku [1], I. I. Orabi and G. Ahemadi [59]. Hysteretic type nonlinear damping 
has been studied by Y. K. Wen [70] [71]. An excellent survey is given by S. 
H. Crandall [33]. These works have primarily focused on the study of the 
response of nonlinear damping model to random excitations. Generally, three 
methods are used in the analysis of nonlinear damping systems under random 
excitations: 

1. The Fokker-Planck approach; 

2. The perturbation approach; 

3. The equivalent (stochastic) linearization approach. 

The main advantage of Fokker-Planck method over all the others is that, the- 
oretically, exact solutions may be obtained when the excitations are Gaussian 
white noise. Unfortunately, its use is limited because of the severe restrictions 
that must be placed on the form of nonlinearities and on the spectral density 
matrix of the excitations. For a more detailed study, refer to [25]. 

If the dynamical system has weak nonlinearities, then the approximate random 
response may be obtained using the classical perturbation theory. First devel- 
oped by Crandall [34], the approach has been later generalized to multi-degree 
of freedom systems by Tung [68]. 

Among the methods mentioned above, the equivalent linearization technique 
has the widest range of applicability. Basically, the method is the statistical 
extension of Krylov Bogoliubov approximation method [49]. Although the 
equivalent linearization method is widely used, it is incapable of displaying 
the nature of nonlinearity. 
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• Infinite Dimensional Model 


To the best of my knowledge, research work on nonlinear damping in a dis- 
tributed parameter system is very limited. 

A. V. Balakrishnan [13], for the first time, established energy nonlinear damp- 
ing model for distributed parameter systems. The damping consists of asymp- 
totically proportional linear damping term and energy nonlinear damping 
term. 

A. Lunardi [55] considered the transverse deflection of an extensible beam 
with hinged ends and the nonlinear damping term considered is nonlinear vis- 
cous damping. By reformulating the nonlinear damping model as a semilinear 
abstract parabolic initial value problem, the author studied the stability and 
the instability of all the stationary solutions and of small periodic orbits near 
stationary solutions for the various ranges of the associated parameter in the 
model. 

H. K Wang and G. Chen [69] considered a vibrating string with one end fixed 
and the other end is installed on a nonlinear damping device whose velocity- 
frictional force relationship as determined by material testing is given by 

T dyj0 2 t) _ 

dx dt 

where f(x) is a a multivalued function. The authors used the method of 
characteristics and nonlinear semigroup theory to study the effect of nonlinear 
boundary damping and analyzed the asymptotic behavior of the solutions of 
such systems. The u>-limit set of the dynamical system and the asymptotic 
rates of various solutions to the u>-limit set are determined. 

Linear Damping (Distributed Parameter Systems) 
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• Strictly proportional and asymptotic proportional damping 

To formulate the internal passive damping, strictly proportional and asymp- 
totic proportional damping operators have been reported in A. V. Balakr- 
ishnan [10] [11], as well as in G. Chen and D. Russell [28], S. Chen and R. 
Triggiani [32]. F. Huang [42] [43] studied the spectral property of the systems 
in the form 

x(t ) + Bx(t) + Ax(t) = 0 

where B is a closed linear operator related in various ways to A a with 1/2 < 
q < 1. The author obtained some fundamental results for the holomorphic 
property and the exponential stability of the semigroups associated with these 
systems. 

Strictly proportional damping operator is essentially the square root of the 
stiffness operator A. In this case, the eigenvalues have the proportionality 
property 

*(A n ) ( 

3(An) 

The drawback of strictly proportional damping is that the damping opera- 
tor contains nonlocal feature, which is unnatural if we consider that internal 
passive damping is due to the structure material itself. However, if strict 
proportionality is relaxed to asymptotic proportionality, i.e. 

*(An) . . 

lim — — = constant 

n-oo $(A n ) 

then, the nonlocal feature can be avoided. 

A. V. Balakrishnan, based upon his theory on the fractional power of closed 
linear operators [14], explicitly calculated the strictly proportional and asymp- 
totic proportional damping operators for the beam bending model [11], in 




\ 
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which one end of the beam is clamped and the other end has an end-body 
attached to it. In [10], the strictly proportional damping operator is given 
explicitly for beam torsion model. 

• Boundary Damping 

There has been considerable amount of literature on the distributed parameter 
systems with boundary damping (including boundary active control). G. Chen 
[29] [30]. J. Lagnese [50] [51], and R. Triggiani [52] [67], to name a few, have 
studied the energy decay of solutions of wave equation on bounded domain 
with damping only on the boundary. In general, in this kind of study, the 
concrete question to be answered is: 

Under what conditions is it true that there is an exponential decay rate for 
E(t), i.e. 

E(t ) < Ce-^EiO), t > 0 

for some positive w? In the above, E(t) stands for the total energy of the 
vibrating systems, which needs to be properly defined. 

• Interior Point Damping 

G. Chen, M. Coleman and H. H. West [31] and K. Liu [54] studied the energy 
decay rate of a coupled vibrating string with a point damper installed at the 
coupling point. In these studies, it is assumed that a damper applies damping 
at an interior point where two strings couple. Possible mechanical designs are 
also proposed in [31]. 

As far as linear damping problem in finite dimensional space, certain work has 
been done. Modern computer-based techniques (such as finite element method) 
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enable the structural engineers to make highly precise calculations of the mass and 
stiffness matrices of elastic structures. These calculations, in turn, lead to quite 
accurate estimates of natural frequencies and mode shapes, particularly for the lower 
modes. Structural damping properties, on the other hand, tend to be much less 
accurately calculated; indeed, one usually simply guesses at modal damping factors. 
Limited work has been done in this area [2] [3] [40]. D. F. Golla and P. C. Hughes [40] 
developed a method of constructing linear damping matrix for viscoelastic structures 
in the framework of finite element method. They assume that certain material 
constants are available (ultimately by measurements) for each constituent material 
of the structure, just as known densities determine mass properties, and known 
elastic constants, such as Young’s modulus, determine stiffness properties. These 
measured viscoelastic material constants permit a set of equations to be formulated 
for the dissipative properties over all parts of the structure. The method merges 
naturally with finite element method and is a natural extension of it. 

1.2 Objectives and Contributions of the Disser- 
tation 

Experimental data has clearly indicated the nonlinear nature of the internal friction 
damping of large flexible space structures. Then, what is the frequency response of 
a nonlinearly damped structure? Until now , Monte Carlo simulation has been the 
only method of computing the frequency response, due to the lack of the parallel 
theory as in linear systems, in which the frequency response can simply be obtained 
from transfer functions. This dissertation provides an analytical method of com- 
puting the frequency response of single DOF oscillator with nonlinear damping. In 
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spite of its fund amen tfil importance, the nature of internal damping has been little 
known. In the modeling aspect, this dissertation proposes an energy type nonlinear 
damping model and the corresponding stationary probability density with white 
noise input can be obtained explicitly. Theorem 3 gives an interesting result, in 
terms of Krylov-Bogoliubov approximation, concerning the modeling and identifi- 
cation of nonlinear internal damping in flexible structures. This work also serves a 
contribution to the random vibration theory by providing a method of computing 
the first and the second order statistics (steady state probability density, correlation 
function and spectral density) of nonlinearly damped oscillators with white noise 
input. 

A Krylov-Bogoliubov type approximation is established for systems having in- 
finite number of DOF’s and its error estimate is obtained. Comparisons are made 
between nonlinear damping (linear stiffness) models and nonlinear stiffness (linear 
damping) models, and between nonlinear damping (linear stiffness) models and lin- 
ear models. 

A group of sufficient conditions for strong stabilizability is provided for general 
distributed parameter oscillation system, taking the actuator saturation into con- 
sideration. These are the weakest sufficient conditions obtained so far and it is 
found that the nature of internal damping is not crucial in guaranteeing the strong 
stabilizability. 

By extending the notions of characteristic equation and its root locus to our 
distributed parameter oscillation system, we studied the nature of active damping. 
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1.3 Organizations of the Dissertation 

This dissertation is organized in the order of dimensions of the model: 

1 . Chapters 2 and 3 are on the single degree of freedom nonlinear damping mod- 
els; 

2. Chapters 4 and 5 are on the multi- but finite- degree of freedom nonlinear 
damping models; 

3. Chapters 6 is on infinite dimensional energy type nonlinear passive damping 
models; 

4. Chapter 7 is on the active damping of a uniform Euler-Bernoulli beam with 
one end clamped and the other end free, with a tip mass. The notions of 
Characteristic equation and Root locus are extended to distributed parameter 
systems. 

5. Chapter 8 is a summary of conclusions and a list of some related open prob- 
lems. 
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Chapter 2 

The Spectral Density of 
Nonlinear Damping Model: 
Single DOF Case 


2.1 Introduction 

The problem of characterizing the damping mechanism in flexible structures has 
received renewed attention in recent years in connection with the need to stablize 
flexible flight structures such as antennas deployed in space. Experimental evidence 
suggests the need for nonlinear damping model and the need to consider the effect 
of random disturbances due to the uncertainties in system parameters and the en- 
vironment. One of the most important subjects in nonlinear random vibration is 
to obtain the second order statistics, i.e., correlation function and spectral density 
of the stationary response, because they provide average amplitude and frequency 
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information about the sample histories. Unfortunately, up to now, the only practi- 
cal method available is Monte Carlo simulation and there is no analytical technique 
for the second order statistics of nonlinear systems [33]. This paper presents an 
analytical technique for computing correlation function and spectral density of the 
stationary response of nonlinear damping model subject to white noise excitation. 
The basic nonlinear damping model we consider is 

x + 2£o> 0 x -f 7_D(x,i) +u>oX = an(t) (1) 

where 7 > 0 is a small constant because the damping in flexible space structures, 
whatever its nature, is small. 

The corresponding Fokker- Planck equation is given by 

Tt = - y Tx + + 2f “ w + 

d 

= L oP + 1-r^[D(x,y)p] 

limp(<, x, ylzo, Vo) = 6{x - x 0 )6(y - y 0 ) (2) 

Notations 

z(t)H D(x(t),y( ()); 

Ps( x ,y)' stationary density of (a:(<), y(<)), i.e., the invariant measure; 
m iJ = fl K> x 'y } Ps( x ,y)dxdy , i,j m 0, 1,2, • • •; 
p(t,x,y\x 0 ,y 0 ) : the fundamental solution of (2); 

Po(<,z,y|zo, yo) : the fundamental solution of (2) with 7 = 0; 

T{t) = exp(Lot); 

9(t,s,x,j/|x 0 ,y 0 ) = If niD(u,v)po(t — s,x,y|u,u)p(s, u, u|x 0 , yo)dudv. 
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It is well-known that po(t, x, y|xo, yo) is a two-dimensional Gaussian density func- 
tion. Its mean vector and covariance matrix can be found by straightforward calcu- 
lation, as 


1 »«x(0 ^ 
V 7 »y(0 / 


t At 


) 



\ yo 

J 


e-**‘ f 

U-’o cos(u>„2 — 0) 

sin w„t 

Lu' n 


— u’o sinu>„* 

u> 0 cos(u> n t -f 0) 



( 3 ) 


S(<) = r 


( 0 


\ 


*l>(0 *J(0 ) 




4*1 


0 1 


e -2{u;o« 


sin(2w n i - 0) + 1] 


1 — cos 2u: n t 


1 — cos2u> n i 

^[^sin(2a;„t + 0 ) - 1] 


( 4 ) 


where 


«n = Uoy/l - e 


0 = tan 


i 


Later on. we will need the notation 


So" f lim m = ~ A j- 

t — oo 4£u.’o 


3 


0 1 


and without ambiguity we will often denote E(f) by £. 
Assumptions on Djx^y) 

(Al) D(x,y) is differentiable with respect to y; 


( 5 ) 
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(A2) 3 A' > O.fc > 0 such that \D(x,y)\ < I<[ 1 + (x 2 + y 2 )*] for (x,y) 6 JR 2 ; 
(A3) niij are finite for all nonnegative integers i,j. 

Of course, to satisfy the energy nonincrease requirement, we also need 

2^0 y 2 + D(x, y)y > 0 (x, y) e JR 2 


2.2 Results 

2.2.1 An equation of spectral density 

Lemma 1 Under assumption (A2) and (AS), it holds 

r*+j£.cc ( ^ L s ’ *' y l- To ’ 2/0 ) = 0 v 0 < s < A (x 0 ,yo) € JR 2 (6) 

provided 

x2+y”-ec /, ^ ,a '’ 2/ ^ 0 ’ = ° > °’ ( x o»J/o) € JR 2 (7) 

Proof: First, by Schwarz inequality, we have 

|ry(As,.r,t/|xo,t/o)| < JJ ^JD(u,v)\p 0 (t — s,x,y\u,v)p(s,u,v\x 0 ,y 0 )dudv 

- [ JJ R2 \ D ( U > v)\ 2 Po (t - s , X, y|tt, v)p(s, it, v|x 0 , y 0 )dudv] 1/2 

X [ JJ ~ ’ S ’ X ' y \ u ' v ^ s ' u ’ v l x °’ yo)dudv] 1/2 

The second sqare root term goes to 0 as x 2 + y 2 — ► oo by the assumption (7). 
Therefore it is sufficient to show that the first sqare root term is bounded for all 
(3‘,?/) G //?-’. 

In fact, we have 

JJ ^JD(ii,v\ 2 p 0 (t - s,x,y\u,v)p(s,u,v\x 0 ,y 0 )dudv 
- 2 TJkkft - tJI K 1 Wu ' v '>^ s ’ u ’ v \^,yo)dudv (8) 
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In addition, by assumptions (A2) and (A3), we have 

[f f[ \D{x,y)\ 2 p{t,x,y\x Q ,yo)ps{xo,yo)dx 0 dyodxdy 

JJ mJJ nv 

= [[ \D{x,y)\ 2 p s {x,y)dxdy 

JJ w 

< JJ ^K 2 [l + ( x 2 + y 2 ) k ] 2 p 3 {x,y)dxdy < oo 

By Tontlli ' s lemma, we know that 

\D{x,y)\ 2 p(t,x,y\x 0 ,yo)p>(xo,yo) € L'ilR 2 <8> 1R 2 ) 

Then by Fubinr s theorem, we know that 

/ / |D(a’, t/)| 2 p(t, x, y|-, -)dxdyp s (-, •) € L l (JR 2 ) 

JJ JR* 

which implies |D(-,-)IV(<i'»-|*o,yo) € I 1 ^ 2 )- Therefore, by (8), the first sqare 
root term is bounded for all (x, y) £ HI 2 - ^ 

Theorem 1 Under the assumptions (Al)-(AS), the following equations hold: 


$xx(“>) = 


4£u$m 2 ,o 


+ 2 7 £[ 




u} 2 — u>q + i2£u>ou: 


] 


(u;2 - u> 2 ) 2 + 4 t 2 ufc 2 

a / \ 4^ 0 7no,2^ 2 , o„, .cvr i 

“ (u> 2 - u.- 2 ) 2 + 4£ a wgu> 3 +27 L*-u£ + »2£wbJ 


( 9 ) 

( 10 ) 


U'htr( 


'i’x ;{u>) = [ Rxz( t )e twr dr 

Jo 

roo 

^y..(w) = / 

v 0 

a?iJ 5ft. Cl denote real, imaginary parts. 

Proof:( 2) is equivalent to the following integral equation 

p(t,r,y\xo,yo) = Po(l,x,y\x 0 ,yo) 

t 

+7 J q T {t - s)-q-[D(x, y)p(s, x, y |x 0 , yo))ds 
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or more explicitly, 


p{t,x,y\x 0 ,y 0 ) = Po(t,x,y\x 0 ,y 0 ) 

+1 J q JJ ^Poit - s,x,y\u,v)-^-[D(u,v)p(s,u,v\x 0 ,y 0 )]dudvds (11) 


Performing integration by part in the second term and noticing the relation 


0} 


e -(uj 0 t 


n (t,x,ylu,v) = 
0 v 


.SPO,. I \ 
sinw n t— (t,x, j/|u,u) 


we obtain the following 


~^===^cos (u n t + 0)-^-(t, y|u, u) 


1 L II ~ s ' x 'y\ u ' v )fo{E > ( u ’ v )P( s ’ u ’ v \ x Oyyo)]dudvds 
1 J 0 JJ ~ ■s,x,y\u,v)D(u,v)p(s,u,v\x 0 ,yo)dudvds 

Q rt e -{u/o(f-s) 

= ”1 777 / smu n (t - s)q(t,s,x,y\x 0 ,y 0 )ds 

(JJ J 0 i jJ n 

Q rt t 

+1 0y Jo cos E^«(f - a) + 0]q(t, s, x, y|x 0 , y 0 )ds 

Thcrelore, we obtain the following integro-differential equation for p(t,x,y\xo, Vo)- 


p(t,x,y\x 0 ,y 0 ) = Po{t,x,y\x 0 ,yo) 

Q rt e -(vo(t-3) 

+7 <9t Jo sinw n(* - s)q(t,s,x,y\x 0 ,y 0 )ds 

Q rt e -tw 0 (<-i) 

+1 'dy Jo cos[u n (f - s) + 0]q(t,s,x,y\x o ,y o )ds(l2) 

Based upon (A2) and hence Lemma 1, we find that the marginal transition 
probability density p[i , .r|.r 0 . t/ 0 ) satisfies the following equation, by integrating (12) 
with respect to y over JR 1 , 

d f l e -^o(<-*) 

p[t,x\xo,yo) = Po{i,z\xo,yo) + 7(-z- sinu; n (f — s) 

ox Jo u n 

X II ~ s ' x \ u ' v ) D ( u ’ v )P( s ’ u ' v \ x <>iyo)dudvds (13) 
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After integration, the third term vanishes by Lemma 1. And in (13), po{t, ®|xo> I/o) 
is a one-dimensional Gaussian density with mean 

[u^o cos(u>„< — 6) + y 0 sinu; n <] 

UJ n 


and variance given by 

2 ~—2(wot 

17-sll- f— sin(2w n ( - 0) + 1)] 

Multiplying (13) by xoxp s (xo,yo) and integrating with respect to (s,£o, yo) over 
Hi 3 , we obtain, for t > 0, 


Rrr(t) = [[ x 0 [ xpo{t,x\xo,yo)dxp s (xo,yo)dx 0 dyo 

JJ B? J - oo 

-1 /„' 8in “'" ( ‘ " s) ILJJ*>L- Mt ~ 


xD(u, v)p(s, u, r|x 0 , 2/o)zoPs(zo, yo)dudvdx 0 dy 0 ds 


= // x 0 - [u.' 0 x 0 cos (u n t - 0) + yosinu; n t]p,(xo, yo)dx 0 dy 0 

JJ n? u?„ 

rt 

-7 / sinu? n (< - s) 

Jo U n 

x // ff x 0 D(x,y)p{s,x,y\xo,yo)p 3 {xo,yo)dx 0 dyodxdyds 
JJ B? JJ B? 




= ™2,0 


= cos(u : n i - 0) 


VT^i 1 

rt 

—7 / silUJ n (< — 5)/2 X2 (s )£/6 

7o U.'n 


(14) 


In the above, we used the fact that 


m 


1 , 1 = // xoyoPs{xo,yo)dxody Q = 0 
JJ B? 


which is generally true in nonlinear mechanics due to symmetry. If m^i ^ 0, we can 
get similar result with similar proof, and the only difference is an extra term with 
mip as coefficient on the right of (9) and (10). 
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Since R xx (r) is even, hence 

/cc roo 

R xx (T)e' UT dT = 2 R xx (t) coswrdr (15) 

And in general, we have the following relation 

/ / f(t — s)g(s)ds cos utdt 

Jo Jo 

f ,>z ‘ roc roc 

= / f{i)cosu:tdt g(t) cosu>tdt — I f(t)s\nutdt g(t) sinutdt (16) 

♦'O Jo Jo / 

if each of the integrals in (16) exists. Also, from the identity 


i 


CO — 0* 




sin u: n te xu,t dt = 


(17) 


— u> 2 — i2£u>u Q 

one can easily obtain two identities corresponding to the real and imaginary parts 
of both sides. 

Next, performing Fourier transform on both sides of (14) using (15), (16) and 
(17), one obtains 


$**(w) = 


4{u?§m 2 ,o 


(a,’ 2 - u.q ) 2 + A£ 2 u:lu 2 


I (^ ~ Rxz{t) COS uitdt 4- 2£uqlj / 0 °° R xz (t) sin u>tdt 

(w 2 — u?q) 2 + 4£ 2 u$u 2 

4gu’gm a,o 0 »[(^ 2 - ul - ;2£w 0 u>)tf«(u>)] 

(w 2 - w 2 ) 2 + 4£ 2 w 2 u> 2 + 7 (o> 2 -u;g) 2 + 4£ 2 u;gw 2 


4£u.’g?7? 2 ,o 


+ 27*[- 


*«(w) 


( u '' 2 — ^'o) 2 + 4£ 2 o,’quj 2 ' L u 2 — u>o 4- i2£utou? 

For the second identity (10), one can similarly first show 

e -Zu 0 t 

Ryy{t) = 77?o.2-7==COs(u; n f + 0) 

vi-r 

ft e -{wo(*-») 


-7 


s: 




cos[a> n (t - a) + 0)Ry Z (s)ds 


from which (10) follows. 


(18) 

□ 
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In linear case, i.e. 7 = 0 in (1), it is well known that the spectral densities of x 

and x have the following relation 


U 2 $ XI (U) = $yy{u>) 


The natural question to ask is whether we have similar relation for nonlinear damp- 
ing model (1)- The answer to this question is provided in the following 

Corollary 1 Under the assumptions (Al)-(AS), it holds for (1) 


w 2 $xxM = 


31 


, Mw) 


UJ 2 — U/Q T z2^u.’ol*/’ 


] = 




u> 2 — u>q + i2£u>ou> 


u,' 0 m 2i o = ^ 0,2 


( 19 ) 

( 20 ) 
( 21 ) 


Proof: (19) follows from 


For (20), first one has 


R yy {r) = — Ey(t)x{t + r) 


c P 


dr 2 

d 2 

dr 2 


Ex(t — r)x(t) 
Rxx{t) 


-r-Rxzir) = -Ey(t - r)z(t) 
dr 

= -Ryz(r) 


Then integration by part gives 


»xs(w) = 


_j_ r dRx ^ii e ^dT 

iu> Jo dr 
= - T Ry^e^dr 

ZD JO 


= -1VM 

ID 


( 22 ) 
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in which we have assumed that 


Rx:(0) = JJ ^xD(x,y)p t (x,y)dxdy = 0 
lim R xz (t) = 0 

Therefore, from (22), we have 

<qt ^^(^0 i __ cvr iu>ty xx (u>) . 

u 2 — u>l 4 z2fu> 0 u> u; 2 — 4 i2(u>ou> 

#«(u>) 


= u%[ 


u> 2 — l>q 4 i2{u qu 


from which (10) becomes 


<b (u:) = 46,'q 777 0 , 2 U- 2 „ 2yf **»M 

y ' J (w 2 — u.'^) 2 4 4f 2 u^w 2 “ u> 2 — u>o + 

4£u*m°.*;* (w) 4 ^o 3 ™ 2 .o l 

' 1 xar V ) / o , .2\2 i >1^9. .2 .9 j 


(U.’ 2 - u.^) 2 4 A^USqU 2 


{u 2 - u$) 2 4 4f 2 u>3u; 2 ' 


= - 


4(fa>o(^’o m 2,o - mo, 2 V 2 


(u 2 - w 2 ) 2 4 4£ 2 u, 2 u> 2 
Comparing (23) with (19), one can immediately obtains (21). 


(23) 


□ 


2.2.2 The spectral density 

In this subsection, we establish the perturbation formula of the spectral density of 
(1) with 0(7 2 ) accuracy. First, we need two lemmas. 

Lemma 2 The following matrix relations hold 
£(<) = S 0 - e At Soe ATt 

S' 1 - E -1 e j4l (Eo 1 4 e ATf 'L~ 1 e At )~ 1 e ATt %~ 1 = Eq 1 
(So 1 4 e' ,Tt E- 1 e 4t )- 1 e AT ‘ = Eoe^Ej^E 
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(24) 

(25) 

(26) 



Proof: Define 


E (t) = f E 0 - ^‘Soe 


Atrr^ A T t 


To show (24) is equivalent to showing E(t) = E(t). 
In fact, since 


*«>-/. 


( 


.*(<-») 


0 0 


z A T (t-s) ds 


v° 


we know that E(f) satisfies the following linear differential equation 


2s(l) = >lE(() + S(f)>l T + 

at 


( o o'' 

^0 O 2 ) 


And. based upon the simple fact 


-AEn - E 0 A T = 


0 0 




\vc also have 


(27) 


4S(<) = /1(E(() - Eo) + (£(0- 

at 


= At(t) + t(t)A T + 


0 0 

0 <r 2 


(28) 


Obviously. E(<) and E (/) have the same initial condition 

E(0) = E(0) = 0 

Therefore, the uniqueness of the solutions of linear differential equations implies 

E(t) = E (t). 

By the matrix inversion formula 

(71A + M 2 M 3 Mt)~' = Jl/f 1 - 1 (29) 
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we have, noticing (24), 


E" 1 - E-V^Eq 1 + e ATt H~ 1 e At )~ 1 e ATt E ~ 1 
= (E + e^Eoe ^*)- 1 

— V- 1 

— ^0 

For the proof of (26), we again use (29) and (24) to obtain 

(Ep 1 + e ATt 'E~ 1 e At )~ 1 e ATt 
= [E 0 - E 0 c^ Tl (E + e^Eoe^-V'Eoje^ 
= Eoe^CZ-Eo'e^Eoe^) 

= Eoe^fJ-Eo^Eo-E)] 

= Eo^'E^E 


□ 


Lemma 3 For linear damping model 


x 4- 2£u> 0 i + u)qX = <rn(t ) 


we have 


< x(t) ' 

^ y(t) j 


D(x(i + T ),y(t + T)) = Eoe^Eo 


^rn-l ( 


V’v 


(30) 


where 


4>x 

\ V'v 


= H«, D ^y) ( * ) + y 2 )]dxdy 


( 31 ) 


Proof: Let us use the following notations: 


X = I I *0 = 

y I \ yo 


def I X 0 
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By (25) and (26), we have 

(A' - e A, Xo) T t-\X - e M X 0 ) + A 0 T E^A 0 
= A r [E-‘ - E-V(E 0 -' + c A 7 , E- , e /u )-'e AT '£- l IX 

+[A„ - (e AT, E-'e AI + S;')-'e ATt E-'X] T (e AT <E-'c M + EJ 1 ) 
x[Ao - ( e AT, Z-'e M + Ej 1 )- W 1 *) 

= X t Eq 1 X 

+(X 0 - Eoe^'Eo 1 X) T {e ATt E~ l e At + Zo'){X 0 - Zoe^Z^Z 


Then 


E f D(i(i + r),y(i+T)) 

V y(t) ) 

= // [[ X 0 D(x,y) i =exp[— 1/2(X — e Ar Xq) t Z~^ {t){X — 

JJ wJJ f? 27t v /|E(t)| 

x i=exp(-l/2^ 0 7 ’Eo 1 Xo)d|X|d|Xo| 

2ny/\Zo\ 

= IL D(x ' v) ^^ ex?( - ll2XTzz ' x) SI ^ x ° 

x } exp[— l/2(Ap — Ept A T So 1 X) r (Eo 1 + e A T E 1 e j4T ) 

x(X 0 -E 0 e ATr E^X)]d\X 0 \d\X\ 


= [I D(x,y) L=exp(-1/2A T S 0 ->A) 

JJ * 2 2 tt v /| 5 :o| 

Eqc aTt Eq 1 X 


dxdy 


v /|E(r)||E^ 1 +e' ir '£->(T)^'| 

= E„e' ,r 'E 0 -> // XD(x,y) !—exp(-l/2X T £;'X )dxdy 

JJ R 1 2tt^|Eo| 


e AT X 0 )] 
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In the above, we used the fact that 


lE^HEo 1 + c^ T ‘E- 1 (0e >lt | = 1 


which could be easily verified by using (26). □ 

Theorem 2 Under assumptions (Al)-(AS), the spectral density of (1) is given by 

4£wgm 2i0 


= 


(a; 2 - a; 2 ) 2 + 4£ 2 u> 2 u; 2 


0 . + 4£ 2 u, 2 u, 2 ) + 'I’vK" 2 ~ *1? ~ 4£ 

+ 7 P? - ^o) 2 + 4£ 2 u>ou; 2 ] 2 + 0(7 } (32) 


Proof: If we write 


R*.{r) 

Ry.(r) 


«&(r) W B§)(t) 

«(’■) J y I BgKr) 


+ 


then (Bfffr), ^(r)) is given by (30), i.e., 


*£'M 

C(r) 


= Eoe 




0* 

Ipy 




( OJQ cos(a> n t - 6) - sin u n t \ / if> x 

^ u;£sina ) n t u> 0 cos(u} n t + 6) J \ rp y 


Through one-sided Fourier transform, we obtain 

*(0) fw x = ~ 2 ^o) ± tv 

” K ’ u 2 -ul + i^wow 


(33) 


( 34 ) 


which, together with (9), gives (32). 



EXAMPLE 1: Consider the nonlinear damping model with 


£>(x,y) = i 2m |ar|V n+1 |y| 0 


where m, n are nonnegative integers and 0 < a, 0 < 1. 

Since Z)(x,y) is an even function of x, V’i = 0- And computation gives 


* = ^<2^>” + ” +,+ ^ r(m + i T 5)r(n + 1 + ^ 


1 + 0 , 






EXAMPLE 2: Consider the following saturation type active damping model, 

x + 2 tu> 0 x + 7 tan -1 bx + wjjx = crn(t ) 


(35) 


Similarly, we have rp x = 0. And 

r°° . _i , v^o 


2 £o>o 


J-oo yJ'KO a 


ob 


T b f°° 1 , 

7r£w 0 -/-oo 1 + 6 2 y 2 6XP 

- M - )le^ 

2- v /2£w 0 bo 




where the last equality is based upon the identity [40, p931] 


2x , j r°° e~ t2yi 

*(* y ) = 1 - T e ‘ ' I ¥T^ dt x ' y -° 


(36) 


in which 


*(x) d ^ -7= f X e-* 7 dt 

y/T JO 

It is not surprising to see that the second order statistics depends upon the first 
order statistics m.j. Therefore it is important to compute those numbers m.j, which 
is our next topic. 
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Chapter 3 


On the Stationary Probability 
Density: Single DOF Case 


3.1 The Stationary Probability Density in a Par- 
ticular Case 

To determine rrijj, it is desirable to have p,(r, y) in explicit form. However, to obtain 
p s (x,y) explicitly is itself a very difficult task which has attracted large amount of 
literature. Many researchers have tried to clarify the largest class of nonlinear 
damping models for which an explicit stationary density can be obtained. The most 
recent and most inclusive account of this subject can be found in [22] and [23], 
which include parametric excitation as well as external excitation cases. However, 
the only useful result so far for the linear stiffness, nonlinear damping model with 
only external excitations is the following [37]: if the nonlinear damping model is of 
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the form 


- , ,w£r 2 + x 2 N . 2 


)x + OJqX = <7 Tl(t) 


then the stationary first order density is given by 

P*{x,y ) = Cexp[ J fi(z)dz ] 

^ = ~ J 0 Jo ^ dz ^ d P 

O u) o Jo Jo 


Next, let us consider the case D(x,y) = )y. In this case, the stationary 

density of (1) is given by 

p,(x, y) = C eX p[-^ • ~ y - - J q ^ ^ p(*)<fe] (37) 

From (37) one can realize the following: 

1. if /i(-) is a polynomial such that p(x) > 0 for x large, then all the assumptions 
(A1)-(A3) are automatically satisfied; 

2. after reaching stationarity, x and x are uncorrelated, but are not independent 
because p(x,x) is not separable unless in the trivial case fi{E) = const., which 
reduces to a linear damping problem. 

For later convenience, let us introduce the notation 

m(k) = f f x*exp [— — j^x — ^ / n(z)dz\dx 
Jo a 2 a 2 Jo 

& = 0 , 1 , 2 , * " 

Then the normalizing constant C can be written as 

1/C = — m(0) 

Uq 
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Through substitution of variables 


{ x — l/o> o pcos0 
y = psin0 

we can write 


m «,j 

C f 2 * 


= f cos’ 0 sin-' 6d8 f {2x) i ^ L exp[— ^ f fi(z)dz]d 

Uln JO Jo O V Jo 


U 0 

f 


= < 


bothiandy 


either i or j is odd 

are even 


(38) 


In (38), we have used the following identity 


f 2n ■ 

/ cos’ 6 sin J 

Jo 


ede 


either i or j is odd 


2B(^-, both i and j are even 

where is a Beta function which has the following relation with Gamma func- 

tion 

B(x, y) • r(i)r( y )/r(i + si) 

from which (38) follows. 

In addition, (38) immediately implies that 


u>£m, j = i,j = 0,1,2, 


from which one can be obtained from the other. 

Therefore, the problem of computing m.j becomes that of evaluating the inte- 
grals m(0) and those necessary m(^) for t, j both even. 
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Notice that (32) is not exactly a perturbation expression because the first term 
contains 77*2,0 which depends on 7. In this particular case that D(x,y ) being of 
the form /i( u ' pI 2 ' l ' v )y, we can obtain the perturbation expression of 77*2,0 and hence 
obtain the perturbation expression of 


Corollary 2 Suppose D(x,y) is of the form /x( u '° r , +v )y, then 


m 2,o — 


$xxM = 


<y 2 - 2yipy 


+ 0 ( 7 2 ) 


^(A^uquj ) 2 


-7 




+ 0 ( 7 2 ) 


Proof: First, computation gives 




where 


/, = [ x : [ p(z)dz exp(- 
J 0 Jo 


4fw 0 


x)dx j — 0, 1,2, 


Successive integration by part gives 


(39) 


(40) 


(41) 


'* - W‘ * !/ ° 


t- 2 eoo ^-2 


^ r [ g ( 4k ), s7 (ii)i ' ,(i)exp( ‘^ i)rfi <42> 


4^0 


Therefore, we have 


m(fc) 

m(0) V 4£u;o 


= (•£■)*«-§ 

fc-1 _2 


* r i P-wf^ xtMx)eM -^ x)dx + ° (t!) 
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In particular, 

m (l) 

m( 0 ) 


27 


4£u > 0 4£u> 0 

.4^0 r°°__ ,_ s , 4^u; 0 


X ~T f X H{ X ) QX V( Y-x)dx-\-0{ 7 2 ) 

a 1 Jo o 

a 2 — 2rfip y 


A£u;o 

where we have used the relation 

4£o> 0 


+ 0(l 2 ) 


f xp(x)exp(-^^x)d. 
Jo a* 


which will be proved in Theorem 4- 

Then, by (38), we obtain (39). (40) is easily seen from (39) and (32). □ 

In the case D(x,y ) is not of the form )y, approximation approach has 

to be made. 


\ 


3.2 The Stationary Probability Density in Gen- 
eral Cases 

For simplicity of notation, we surpress the damping term 2£u>oi -f- 7 D(x, x) to, say, 
simply Z)(x,x). Then the nonlinear damping model (1) is rewritten as 

x + D(x, x) + u>lx = <xn(<) (43) 

What we propose is to approximate the exact model (43) by the following mod- 
ified model 

2 2 2 

X -I- n( — 2 * - ~ ) j + = <rn(i) (44) 

where fi(E) minimizes 

/ [D(^^ smtJ),y/2E cosrj?) — fi(E)>/2E cosrJ>] 2 dxl> 

Jo U>Q 
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and is given by 


A r* fOT? 

p(E) = — 7 = / 2 D{ sin xp,y/2E cos xp) cos xpdxp 

7t\/2 E Jo ujq 


(45) 


Presently, the problem of giving a precise estimate of the error of the above 
proposed approximation is still open. But we do have the following side evidences. 

Theorem 3 If D(x,y) is even with respect to x and odd with respect to y, then both 
the exact model (43) and the modified model (44) have the same Krylov-Bogoliubov 
approximation [48]. 

Proof: For the exact model (43), 

daft) 1 ^ 2 ^ 

— ■ — = / D(a sin xp, au 0 cos xp) cos xpdxp (46) 

dt 2nu>o Jo 

By the assumption on D(x,y), 

dm 


dt 


1 f 2* 

/ D(a sin xp,au>o cos xp) sin xpdxp 

2tu>o a Jo 

1 /t 

— I D(a cosxp,auosin ip) cos xpdxp 
Tru> 0 a J-% 


= 0 

For the modified model (44), the Krylov-Bogoliubov approximation is given by 
da(t) 1 f 2 * ,ula 2 


dt 27tu;i 

a ( u o q2 \ 
= -jM— ) 

1 


- f fi(^-)uJoa cos 2 xpdxp 
'o Jo 2 


2iru>o 


J r* 

' D( a sin ip, L>oa cos xp) cos xpdxp 

o 

dm i r 2 * . ,j, 

= / u( — - — )u>oa cos xpsm xpdxp 

2nu> 0 a Jo 2 


dt 2nu 0 a 

= 0 


Therefore, (43) and (44) have the same Krylov-Bogoliubov approximation. □ 
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Theorem 4 If D(x,y) is even with respect to x and odd with respect to y, and let 
V’x, V’y correspond to the exact model (43), V’y correspond to the modified model 
(44) > then 


i>x = = 0 


4£o> 0 [°° , s , 

= ^y = —T xp(x)exp( — -2-x)da 

v G l J 0 <7 


Consequently, 


( x(t) \ ,u>?lx 2 (t + r) + y 2 (t + t) . 

D(x(t + r),y(t + t)) = £ I p{— )y{t + t) 

V v(t) j 2 


Proof: (47) is obvious because both D(x,t/) and the corresponding fi( 0 2 +*-)y are 
even with respect to x. 

For (48), we have the following 


= J J ex p(--^ r2 ) rMdr 

n2 7T Y /2tt y _ * 

— / D( — smip,rcosrp)cosipdrl>rsui 6 
Trio u)q 


2<fu>o / 2£u>o 2\ jflj 

x exp( — r l )rdvar 

<irsT* /T* 


7T (J 

r oo r2ir 


roo r2ir j- 

= / / /)( — sin cos xp)r cos rpdip 

JO Jo UJq 

2£t*>o , 2£u>o 2\ Jftj 

x exp( — r )rdQdr 

-kg* a 1 

= If D(x,y)y^^exp[-^Y-{u>lx 2 + y 2 )]dxdy 

JJ J?2 7T<7* (7* 

= V» v 


The second equality in (48) is easily seen from (50). 
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The solution of the stationary Fokker-Planck equation corresponding to (44) is 
given by 

2 

p 3 (x,y) = Cexp[ 

<J 

where 

1 r 2 fP 

— = — exp[-— I p(z)dz]dp (53) 

c ljq Jo (t Jo 

Example 3: Again, consider the saturation type active damping model (35). 

By the identity 


£ 


/i(z)<iz] 


(52) 



one can compute 


tan 1 (b cos x) cos xdx 


ir VTTv — 1 
2 b 


be 2R 1 


»(E) 

4 


= — -==. / 7 [2^woV^2Ecos + A tan 1 ( b\PlE cos VO] cos 
7rv2E Jo 


= 2£ojq + A 


VT + We - 1 


And consequently, one has 


bE 


-~J 0 ^ z '> dz 

- JT 6 ' /r +2^ 


4_\ 

H — — lnfl + Vl + 26 s E] -f const. 
a 2 b 

The solution of the stationary Fokker-Planck equation: 


p(x,y) = C[ 1 + y/l + 6 2 (u?^x 2 + y 2 )]^ 
xexp[-^^(u;^x 2 + y 2 ) - + fc^x 2 + y 2 )] 


(54) 
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where 


_1_ 

C 


7T ,2fu>(K 

Z^ exp( M ] 

*i (1 + \Zx)^exp[-^^x - ^s/z\dx 


It is easy to realize that p(x,y) achieves maximum at the origin. 

(48) provides a simpler means of computing ip v . From what follows, one can see 
can be easily obtained without employing the integral identity (36) as in Example 

2 . 

By (48) and(54), we have 


i>v = 


4£w 0 f°° Vi + 26 2 x — 1 


= -/ 
b 


■L 


exp(- 


4£u>o 


x)dx 


b * ' &* 

y( y — \)e~ y2 dy exp(^|^) 

' * kv 6 2 <7 2 ' 


b 2 a 2 


(substitution of variable 1 + 2 b 2 x = ^ y ) 


e~ y2 dy exp(^£) (integration by part) 


y/Z&kJ&P 

sfr ° . v/5^ . . 2£ojq N 


which is the same as obtained in Example 2. 
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Chapter 4 


The Spectral Density of 
Nonlinear Damping Model: 
Multi-DOF Case 

4.1 Introduction 

As indicated in the first part of this work [70], the problem of characterizing the 
damping mechanism in flexible structures has received renewed attention in recent 
years in connection with the need to stablize flexible flight structures such as an- 
tennas deployed in space. Experimental evidence [8] suggests the need for nonlinear 
damping model and the need to consider the effect of random disturbances due to the 
uncertainties in system parameters and the environment. One of the most impor- 
tant subjects in nonlinear random vibration is to obtain the second order statistics, 
i.e., correlation function and spectral density of the stationary response, because 
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they provide average amplitude and frequency information about the sample histo- 
ries. Unfortunately, up to now, the only practical method available is Monte Carlo 
simulation and there is no analytical technique for the second order statistics of 
nonlinear systems [32], This paper presents an analytical technique for computing 
correlation function and spectral density of the stationary response of n— DOF non- 
linear damping model subject to white noise excitation. The single DOF case was 
fully investigated in the first part of this work [70]. 

The basic nonlinear damping model we consider is, for x(<) € ZR", 

Mx+ DqX + 7l)(x, x) + Kx = <rn(f) (55) 

where 7 > 0 is a small constant because the damping in flexible space structures, 
whatever its nature, is small. Here, M > 0, D 0 > 0, and n(t) is m dimensional white 
noise, a is an n x m matrix. Since we are only interested in oscillation problem, we 
assume that (55) has no rigid body mode, i.e. K > 0. 

The corresponding Fokker-Planck equation is easily seen to be 

= -y T V x p + V v • [(M^Kx + M~ 1 D 0 y)p] 

+7V„ • ( M~ 1 D(x,y)p ) -f ^tr (M~ l o<r T M~ x V*p) 

= Lop + 7V,, • (M~ 1 D(x,y)p) 

limp(t,x, J /|x 0 ,i/o) = 6(x - x 0 )6(y - y 0 ) (56) 

flotations 

*(<)"£(*(<).»(«)); 

p t (x,y): stationary density of (x(t),y(<)), i.e., the invariant measure; 
p(t,x,y\x 0 ,y 0 ) : the fundamental solution of (2); 
p 0 (t,x,p|x 0 ,y 0 ) : the fundamental solution of (2) with 7 = 0; 
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T(<) = exp(L 0 t); 

q(t, s, x, y|x 0 , y 0 ) = f jfjfj?2nl>( u > v)p 0 (t - s, x,y\u, v)p(s , u, u|x 0 , y 0 )«i|ti |c?|v|. 


It is well-known that p 0 (t,x,y\x 0 , y 0 ) is a 2n-dimensional Gaussian density func- 
tion. Its mean vector and covariance matrix are given by 


/ 


m x (t) 


\ ™v(*) 



( 


(0 *12(0 
^ *21 (0 *22(0 



where the matrix A is defined by 

( 0 


A 

A = 


—M~ l D 0 


2nx2n 


and 


E(*) = f [* e A{t — ] 
Jo 


/ 


0 0 

0 M~ l cra T M -1 
Later on, we will need the notation 


e A T (i-')d £ 




Eoo = f lim £(*) = 

t— mx> 


11 ^12 


Pi Pi 


2n x 2n 


(57) 


(58) 


(59) 


12 r 22 

and without ambiguity we will often denote E (t) by E. 

Assumptions on Z)(x,y) 

(Al) Each component of D(x,y) is differentiable with respect to y; 

(A2) 3 K > 0,k > 0 such that ||L>(x,y)|| < K[ 1 + (||x|| 2 + ||y|| 2 )*] for (x,y) € 2R 2n ; 
(A3) J fpn ( 1 1 x 1 1 2 + ||y|| 2 ) fc p.(*>y)<*M rf lyl < 00 for all nonnegative integers k. 

Of course, to satisfy the energy nonincrease requirement, we also need 


y T D 0 y + 7 y T D(x, y) > 0 (x, y) € lR 2n 
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4.2 An Equation of Spectral Density 

Lemma 4 Under assumptions ( AS ') and (A3), it holds 

.. ,bJ!mb 9(Us,x,y\x 0 ,y 0 ) = 0 VO < s < <, (x 0 ,y 0 ) € IR 2n (60) 

provided 

lim p{t,x,y\xo,yo) = o VOO, (x 0 ,y 0 ) € lR 2n (61) 

PROOF: First, by Schwarz inequality, we have 

||g(<,s,x,y|x 0 ,y 0 )|| 

- / / J?2n ll-^( u ’ t, )IIPo( < - 5,x,y|u,t>)p(s,u,u|i 0 ,y 0 )rfu^ 

“ (/ j ~ s ,x,y\u,v)p(s,u,v\x 0 ,yo)dudv ] 1/2 
X ^ ~ ' S ’ x ' y \ u ' v )p( s ’ u > v l x °’ yo)dudv ] 1/2 

The second square root term goes to 0 as ||x|| 2 + ||y|| 2 — * oo by the assumption 
(61). Therefore it is sufficient to show that the first square root term is bounded for 
all (x, y) € 2R 2n . 

In fact, we have 

/ / /p Jl jC> ( u,t, )l l 2 P°( < ~ «,a:,y|«,v)p(s,u,t;|xo,yo)rf«^v 

In addition, by assumptions (A2) and (A3), we have 

U *?*! . l/l x o, yo)p,(xo, y 0 )dx 0 dy 0 dxdy 
= JJ R 2 n \\D{x,y)\\ 2 p l {x,y)dxdy 

< J j^K 2 [ 1 + (||x|| 2 + ||y|| 2 )fc] 2 p*(x,y)dxdy 

< oo 
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By Tonelli's lemma, we know that 


\\D(x, y)\\ 2 p{t, x, y bo, y 0 )p,(x 0 , y 0 ) € L'iJR 2 ” <g> R 2n ) 

Then, by Fubinrs theorem, we know that 

// \\D{x,y)\\ 2 p{t,x,y\-,‘)dxdyp s {-,-) € L\R 2n ) 

J J #? 2r 

which implies 

l|£(-»‘)ll 2 P(<»vl*0ilfo) € L 1 (R 2n ) 

Therefore, by (62), the first square root term is bounded for all (x,y) € R 2n . Thus 

the claim is proved. ^ 

In the following theorem, we establish an equation for the correlation matrices 

which plays a fundamental role for later development. 

Theorem 5 Under the assumptions (Al)-(AS), the following equation holds: 

R ix (t) = iU(0)*f 1 (T) + ^(0)^(7-) 

-7 f - s)ds (63) 

Jo 

for t > 0 . 

PROOF: (56) is equivalent to the following integral equation 


p(t,x,y\x 0 ,yo) = Po{t,x,y\x 0 ,y 0 ) 

+7 f T(t - <s)V„ • [M~ 1 D(x,y)p(s,x,y\x 0 ,yo)]ds 
Jo 


or more explicitly, 


p(t,x,y\x 0 ,yo) = po{i,z,y\xo,yo) + 1/ J Q J K:in Po{t - s,x,y\u,v) 

xV v [M _1 D(u, u)p(s, u, u|x 0 , yo)]d|u|d|u|ds (64) 
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Performing integration by part in the second term and noticing the relation 


V v p 0 {t - s,x,y|u,t>) = -$f2(*- s ) v *Po(<-s,x,y|u,t>) 

- s)V v p 0 (< - s,x,y|u,v) 


we obtain the following 

/ p 0 (t - s, x, y|tt, v)V v • [M~ l D(u, v)p(s, u, t>|x 0 , y 0 )]d|u|d|v| 

Jfpn 

= - / [V v p 0 (t - s,x,y|u,u)] T Af _1 I>(u,u)p(s,u,v|xo,yo)^|«Mk| 

JR 2n 

= Jjy*# -«s,x,y|u,v)] T $ 1 2(f - s)M~ 1 D(u, v)p(s, u, v|x 0 , y 0 )d|u|d|v| 

+ L„ [^vPo(< - 5,x,y|u,t;)] T $22(< - s)M~ 1 D(u, v)p(s,u, v|x 0 > yo)d|u|d|t>| 

J R~ n 

= V* • / $i 2 (< - s)M~ 1 D(u,v)p 0 (t -s,x,y|u,v)p(s,u,v|x 0 ,y 0 )d|«|d|u| 

JR 7n 

+V v • J $ 22 {t - s)M~ 1 D(u,v)p 0 (t - s, x, y|u, v)p(s, u, u|x 0 , y 0 )d|tt|d|n| 

= V r • [$i 2 (t - s)M~'q(t, s, x, y|x 0 , y 0 ) 

+ V y • [$22(< - s)M~ 1 q(t, s, X, y|x 0 , y 0 ) 

Therefore, we obtain the following integro-differential equation for p(f,x,y|x 0 ,yo): 

p(t,x,y|x 0 ,y 0 ) = Po(<,x,y|x 0 ,y 0 ) 

+7V X - f $i 2 (< - s)M~' l q(t,s,x,y\xo,yo)ds 
Jo 

+7^v • / $22(< - s)Af -1 g(<,s,x,y|x 0 ,y 0 )ds ( 65 ) 

JO 

Based upon (A2) and hence Lemma 1 , we find that the marginal transition 
probability density p(t,x|x 0 , yo) satisfies the following equation, by integrating ( 65 ) 
with respect to y over 2R", 

p(f,x|x 0 ,yo) = Po(f,x|x 0 ,y 0 ) + 7V* • / $ 12 (< - s)M ~ 1 

Jo 

x Ik*" ^ u ’ “ s ’ *l u ’ v )p( s > w|*o» yo)^|«|d|u|c?s (66) 
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After integration, the third term vanishes by Lemma 1. And in (66), p 0 {t, x|x 0 , yo) 
is an n-dimensional Gaussian density with mean 


$n(f)x 0 + ^*12(0^0 

Multiplying (66) by x 0 x T p t {x 0 , y 0 ) and integrating with respect to (x,x 0 ,yo) over 
IR 3 ", we obtain, for t > 0, 

R xx {t) = / x 0 x T Po(<,x|xo,yo)p»(xo,yoM|x|d|x 0 |d|yo| 

Jr 3n 

+ 7 / x 0 x T V I • / * 12 (< - s)M~ i [ D{u,v) Po (t-s,x\u,v) 
jR3n JO JfP n 

xp(s, u, v|x 0 , y 0 )d|u|d|u|<fsp 4 (x 0 , yo)d|x|d|x 0 |d|yo| 

= [ xo[^n(0xo + $i2(<)yo] T p*(xo,yo)rf|x 0 |d|yo| 

J F? n 

-7/ x 0 r / D T (u,v)M- 1 $l 2 (t-s)p 0 {t-s,x\u,v) 

JrP n Jo Jrt 7n 

xp(s, u, v|x 0 , y 0 )<f|u|d|u|p s (x 0 , yo)rf|x|d|x 0 |d|yo!ds 
= ^,(0)^(0 + ^(0)#»(0 

- 7 /7 / xo^tt,^" 1 *^-*) 

Jo JR 7n JR 7n 

xp(s, u, u|x 0 , yo)Pa(xo, yo)d\u\d\v\d\x 0 \d\y 0 \ds 
= J?«r(0)*n(0 + R XV (0)^ 2 (t) 

-7 [* R X2 (s)M-'$* 2 {t - s)ds (67) 

Jo 


The significance of this theorem is that R xx (t ) is expressed in terms of R xz {t), 
while the unknown matrix R XZ { T ) appears only in the first order coefficient of 7* 
Therefore, if one obtains the fcth order perturbation of R X z{ T ), then R xx (t) is im- 
mediately given by (63) with error in the order of 0( 7 fc+2 ). 
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4.3 The Spectral Density 


In this section, we establish the formula of the spectral density of (55) with 0(~f 2 ) 
accuracy. First, we need two lemmas. 


Lemma 5 the pair 


{A, 


0 

M~ l a 


being completely controllable, then the following matrix relations hold 

S(<) « Eoo - e At E oc e ATt 

S' 1 - 4- e AT, Z-'e J> ‘)-'e AT, E-' = E’ 1 

(Sp + c AT, £-'i A, )-'c AT ‘ = E 00 e /,r, E- 1 E 


( 68 ) 

(69) 

(70) 


PROOF: Under the controllability assumption we can know that both E(tf) and E c 
are positive definite for all t > 0. 

Let us define 


2(0 = - e"E»e 

To show (GS) is equivalent to showing E(f) = E (t). 
In fact, since 


i4tv *A T t 


:«) = /' 

Jo 


A (is) 


0 


0 


0 M-'oa T M~ l 


oA T (t-,) di 


we know that E (t) satisfies the following linear differential equation 

0 0 

y 0 M~ l acr T Af -1 


-E(f) = ,4E(f) + E(f)yl T + 


(71) 
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Letting t — ► oc, we obtain the relation 


0 0 \ 
0 M~ x aa T M~ l J 


by which we also have 



A(±{t) - Eoo) + (£(*) - EoM r 


.4E(*) + E(*M r + 



0 ] 

M~ l aa T M 1 J 


Obviously. E(<) and E(2) have the same initial condition 


( 72 ) 


E( 0 ) = E( 0 ) =0 

Therefore, the uniqueness of the solutions of linear differential equations implies 
By the matrix inversion formula 

(Mi + vU 2 .U 3 .U 4)- 1 = M^ - ( 73 ) 


we have, noticing (24), 

£-1 _ v-i e -4<(v-i + e * T 'Yr l e M )- x eA Tty £r l 
= (E + c^'E^e ^)- 1 

_ V-l 
“'eo 

For the proof of (70), we again use (73) and (68) to obtain 

(E- 1 + 

= (Ecc - E 00 e / ' T< (E + e At T, <x> e ATt )~ 1 e At 'L 00 ]e ATt 
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= E.„c aT ‘[1 -Z-!c A, Z„e AT ') 

_ v e ^ T *v-iy 


□ 


Lemma 6 Assume flu matrix pair 




0 )> 

M-'o- j 


is completely controllable. Then for the linear damping model 


M x + Dx + Kx = a n(t) 


wc have 


E 


t 


x(i) 


\ 2/(0 


D r (x(t + T ),y(i + r)) = Eoo^E " 1 


4>x 


2 nxn 


where 


(74) 



PROOF: LcM us use the following notations: 


(75) 


.V d = 
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By (69) and (70), we have 


(A’ - e M X 0 ) T 'Z~\X - t M Ac) + Xl^Xo 
= A^E" 1 - E^e^E" 1 + e AT *E- 1 e ylt )- 1 e >iT ‘E- 1 ]A 
+ [A 0 - (e ATi X-'e At + 
x [A 0 - (e ATi 'Z-'e At + S'^'V^E^A] 

= A' r E-»A 

+ (A 0 - S 00 c ATi ZgX) T {e ATt 2-'e At + E' 1 )^ - 


Then 


x(i) 

E | | P T (a-(t + 7),y(f + r)) 

2/(0 


= [[ f [ X 0 D T (x,y) y 1 -- 

JJ iV” JJ w» 27t v /|E(t)| 

x exp[-l/2(A - t AT Xo) T E~ l (t )(A - e^Ao)] 


x 


277^/1 


L=_exp(-l/2A' 0 T i:' , A'„)rf|A'M|A'„| 




exp[— 1/2(A 0 - E oc e >lTT E- 1 A) T (E; 


x(A 0 - E«.e^ T E- 1 A)]d|A 0 |d|A|} T 


2 r > / iS(rji 
:(A 0 — E 

= [// ^U‘,2/) i=exp(-l/2A r E- 1 A) 

(sr ^A T r^^iy\T ~ 

“ ,<x ’ e ' =<f|x|d|y|] T 


_ V 

— — 'cc 


yiE(r)||Ss> + e» rT E-^(r)e^ j 

^ 1Lx ^ y )^ r { - 1,2XT: 


e A, + £i’) 


E ^' E ' 1 *) 


'T 

0 

? + e j4T ’ r E -1 c j4T ) 


0o 1 A)J|x|</|j/| 
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_ V /» r T v-1 

— ‘-'OC> C ^oo 


Tpx 


V tj v 


In the above, we used the fact that 

IS (OIK + e AT "£,- l {t)t M | = 1 V< > 0 
which could be easily verified by using (70). 

Theorem 6 Under ihe assumption of the complete controllability of 

M.( °, I) 

\ M-'a 

the spectral density matrix of (55) is given by 


¥«(w) - ixH + *; x (u;) 


with 


»rx(w) = 72 a rxr(0)(-D — + ^^(0) Af G(£fa7) 

T 

t /... ' / 




;//crc 


-7 


V o y 


D-iuM -K 1 / G(iw) 
M —iuM /l 0 


■ v-i 


( 


V V\ 


GM + 0( 7 2 ) 


G(iw) = f (-w a Jlf - iwD + if) 


-1 


PROOF: First of all, it should be noted that (63) is valid only for r > 
by the basic property of correlation function matrix 


R,A-r) = RU(t) 


(76) 

° ) 
G(iw) J 

* 0. However, 
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we can easily obtain R x x (r) for r < 0 as 


R xx (r) = R; x (-t) 

Therefore, if we define 

Jo 

then we have the following expression for the spectral density matrix 

*xxM = # xx (u>) + «^(w) 

Next, to find the expression of $„(w), we need to evaluate 

r e At e iut dt 
JO 

To this end. we perform integration on both sides of 

. d 


z iuJt —e At = Ae At e iwt 

dt 


to obtain 


$(?w) = ( — i^l2nx2n A) 


-1 


/ G{iui) 0 \ { D — iuM M 


y 0 G(iu>) 

where the last equality can be easily verified. 
If one writes 

f Rxzi*) 

\ Ryz{s) 


—K —iuM 


/fflw W* 8 ?m | + 

m(>) ) v fijt'w 


then, by Lemma 5, 


R { 2(s) 


\ 


R [ y 0 l(s) ) 


- v e^’E' 1 ' ^ 

— oc e ^oo 




(78) 


( 79 ) 


46 



which implies 


= 

And correspondingly, we have 


/ r \ 


1 nx n 

o 




/ 


i’x 

l V’v 




Z) — iwM 

M 



f G(iu>) 
0 


0 

G?(iw) 



Therefore (77) becomes obvious based upon (63), (78) and (80). 


(80) 

□ 
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Chapter 5 


On the Stationary Probability 
Density: Multi-DOF Case 


5.1 A Necessary and Sufficient Condition of Un- 
correlatedness 

A natural but nontrivial question to ask is in the stationary state, whether x and 
x are uncorrelated as they always are in single DOF case. The answer is no, in 
general. Such an example will be given after the following theorem. Then the next 
immediate question is under what conditions x and x are uncorrelated. Next, we 
concentrate on only linear model 

Mi + Dx + Kx = <rn(t) (81) 

for which we have the following conclusion 

Theorem 7 Assume the matrix A2nx2n is stable . Then, in stationary state , x and 
x are uncorrelated if and only if there exists a self-adjoint, nonnegative matrix P n xn 
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satisfying 


(82) 


KPM = MPK 
DPI< + KPD = aa T 
And in this case, the stationary density is 2n- dimensional Gaussian with mean 
0 and variance 

( P 0 \ 

See = >0 

\ 0 M~ l KP ) 

PROOF: (=>) First we establish the following equation concerning the structure 
of the variance matrix E,*, 


S ao = 


where we use the notation 


P\\ P\2 

, -Pn M-'(KP n -DP l2 ) 


Eoo = 


P\\ P 12 


P\2 P 22 


In fact, by the equation, 

( Jf* 1 4\ Jf\ U\ \ ( 


£ 

dt 


$ll(0 $12(0 

V$2l(0 *»(<) 

we have the following relations 


0 I 

\ -M~'K -M~ l D 


$n (0 $12(0 

$21 (0 $22(0 


612(0 = ^22(0 

6 22 (0 = —M~ l K9i 2 (t) - M -1 Z)$ 22 (0 

Also, by the equation 

/ 

0 

o < j t M~ l 


= /> (° 

J ° \ 0 M~' 

*22(tm £(«) ) 


dt 


(83) 


(84) 

(85) 
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where Q = M -1 <rcr T M -1 , we obtain, in the hght of (84), 

Pi2 = 

JO 

= r *i 2 {t)Q*? 2 {t)dt 

Jo 

= - r^2(t)Q^)dt 

Jo 

— —P T 

— **12 

i.e., P 12 is skrew-symmetric. 

Similarly, by (84) and (85), we have 

Jo 

= r *nit)Q*U*)* 

Jo 

= - r*M(<)c*w(o«tt 

Jo 

= M- 1 KP U + M- l DP? 2 

= M-'iKPn-DPu) 


Noticing that, by (83), 


x and y being uncorrelated 


<=> P\2 — 0 


<=> Soo = 


11 


0 


and that E,* satisfying 


AEqo + Eqo/4 + 


0 


0 M-'KPu 


0 


0 1 


= 0 


then it must be true that 


M~ l KPn = (M^KPuf 

M-'DM-'KPn + (M-'DM-'KPu) 7 = M - V<r T Af -1 


( 86 ) 
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which can be easily reduced to (82). 

(•$=) Suppose P > 0 satisfies (82), then we know 

Eoo = ( P 0 

\ 0 M-'KP 

satisfies (86). On the other hand, the solution of (86) is unique, based upon the 
following general fact [45, p527]: if A mXm and B nXn have no common eigenvalues, 
then the matrix algebraic equation 


^nxm 


A 


m x m 


-B n 


T 

i nXm 


= c n 


has a unique solution T nxm . Therefore, the defined above is indeed the station- 
ary covariance matrix, which is block diagonal. □ 

Now let us consider the following example in which p,(x, y) is not separable, i.e. 
x and y are not uncorrelated. 

Take, for x € -ff? 2 , 


M 

D 

K 


'2x2 


( di 0 

{ 0 d 2 

( “1 0 
{ 0 U , 2 



with d x , d 2 > 0 


with a» 2 ^ 


Obviously the corresponding matrix A is stable. And it is easy to verify that 


{A, 


0 

a 


} 


51 


is completely controllable, i.e. Eqo > 0. 

Next, we observe that for this example, 

KPM = MPK 
=» PK = KP 

*) 

\ 0 P 2 / 

However, on the other hand, 

( d x piu\ 0 \ T 

0 d 2 p 2 ^2 / 

i.e., there is no such ^ 2 x 2 satisfying (82). This, by Theorem #, implies that x and x 
are not uncorrelated in steady state. 

Next, we make the following observations: 


1. (82) is equivalent to 


KPM = MPK 

< 

( -DM~')(MPK ) + ( KPM){-DM~ x ) t + aa T = 0 


(87) 


If D > 0, then -DM~ l is stable because -DM -1 has the same eigenvalues as 
-Af -1/2 Z)M -1/2 which is negative definite. Then the unique solution of (87) 
is given by 


MPK = KPM = r e- DM ~ H ca T e~ M 1 Dt dt (88) 

Jo 

Therefore, the necessary and sufficient condition for x and x to be uncorrelated 
in stationary state is that 


P = M-' 

Jo 

is self-adjoint and nonnegative. 
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2. Let P be the self-adjoint, nonnegative matrix satisfying (82), then P > 0 <=£* 
(DM~',a) is completelly controllable. 

3. In the particular case, D = k 0 cr<T T , where k 0 is a scalar constant, we can 
take P = K~ x in (82). In this case, the stationary probability density is a 
function of the energy E = 1/2 (x T Kx + y T My), 

?.(*.») = (“) n (|V||Al) -1 ^ 3 exp[— + y T My)] 


5.2 Energy Type Nonlinear Damping Model 

In fact, this conclusion can be generalized to the following type nonlinear 
damping model. 


Theorem 8 In the following type nonlinear damping model 


M x p ( 


x T KDx + y T MDy . T 
)o a T Dx -f Kx = (rn{t) 


we assume D is positive definite and commutative with K and M , and (89) is stable. 
Then the stationary probability density is given by 

Ps(x,y) = Ps{E d ) = Cexp(— 2 J ° p(z)dz) (90) 

where 


E d = 1 /2(x r KDx + y r MDy) 

1/C = 2"(\K\\M\\D\Y' /2 ( n 'h) r exp(-2 l%(z)dz)dx (91) 

J =0 Jo •'0 
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with Ij d = fo sin-’ 6dB which has the following iteration relation 


la = 7T 


h = 2 

Ij = j = 2, 3, • • 

PROOF: Assume that p 3 (x,t/) is a function of Ed, then 

V x p,(x,y) = fjJ^ KDx 
V y p,(x,y) = MDy 


Also, we have the relation: 


tr[M"W r M _1 Vjp,(x,y)] = V„ • [A/ W r M 1 V„p,(*,y)] 


Therefore, by noticing the relation 

-y T V I p i + V„ • (M~ l Kxp 3 ) 
= -y T V xPj + x T KM- 1 V yPs 
= 0 


the stationary Fokker-Planck equation 

0 = -y T V lV . + V„ ■ \n(E D )M-'acr T Dyp, + M~'Kx T ,\ + itr(J»f- , <w T M' , Vji>,) 

(92) 


is easily reduced to the following form 

0 = V y • [p{E D )p t M~ x ca T Dy + l/2M"Va T £y^r-] 


or 

0 = V v ■ [M~'a(T T Dy(p{E D )p, + 
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Then it is sufficient for p,(x,y) = p 3 (Ed) to satisfy 

1/2 + h(Ed)Ps = 0 

or, equivalently, 

f^D 

p,(x, y ) = Pt(Eo) = C exp(— 2 / p(z)dz) 

The integrability of p,(x, y) on JR 2n is guaranteed by the fact that both KD and 
MD are positive definite under the assumptions on D. 

The normalizing constant C can be obtained by first making the following vari- 
able transformation 

VKD 0 
0 Vmd 

then changing (u, v) into 2n-dimensional polar coordinate (r,<f> i, • • • , 4> 2 n -i)- O 
REMARK: It is quite natural that in muti-DOF case the stationary probability 
density is a function of Ed instead of E. Consider a simple example, 

m 0 x + Dx + uJqX = a 0 n(t) (93) 




where mo, Wq, <t 0 are positive constants, while D is a positive definite n x n matrix. 

Obviously, the corresponding undamped system is uncoupled and it is easy to 
realize that D > 0 implies (93) is stable. 


In this case, (82) reduces to 

DP A PD = — £/ 

u>& 


therefore, 


and 



p*(x,y) = 


(raou?o ) w/2 

( 7rff o) n 


|D| exp [— 1 /oq (u)qX T Dx + m 0 y T Dy )] 
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5.3 An Illustrative Example 

In this subsection, we consider the following special type of nonlinear damping model 

Mx + (d 0 + 1 n(E))Dx+ Kx = crn(t) (94) 


where D = ac T , d 0 > 0 and ft(E) > 0 for large E. 

We need to find the four n x n matrices ^ x , ifry, Rxx{ 0) a-nd ^?iy(0). 

First, by Theorem 4 we know that the covariance matrix of the corresponding 
linear system is given by 

( I< - 1 0 

Soo = l/2<fo 

l 0 M- 1 


and therefore, by definition, 


t/> v = J ^E)yy T <T<J T QrJmM\eM-do(x T Kx 4 - y 1 
x exp[— d 0 (|l u l| 2 + ||u|| 2 )]d|u|d|u| 


My)]d\x\d\y\ 


lMl 2 + IMl 2 )( 4> r 

2 7T 


= M V cr T f v]n( 

JF? n 

x exp[— <f 0 (||u|| 2 + ||v|| 2 )]d|u|d|v| 

= f dr f d<f> 2 n-i / d<f> 2 n —2 ' ' ' / dfar 2 cos 2 4>\n{r 2 /2) 

Jo Jo Jo Jo 

x (— ) n e _<<or2 r 2n_1 sin 2 ” -2 <f>i ■ • • sin <£ 2 n -2 

7T 

= M'W t [°° x n e~ 2d ° x fi(x)dx 
Jo 7T 

x2/ 0 /l • * ' hn-3{hn-2 ~ hn) 

2(n-l) 


= M 1 ctct t — — TJ Ij 

n * n ,io 
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where 


k„ ^ ( 2d 0 ) n f x n p(x)e 2doX dx 
Jo 

Similarly, we obtain 


V>x = 0 

Rxy{ 0) = 0 

By noticing that the corresponding stationary density p,(x,y) is given by (90) and 
(91), we obtain 

•Rrx(O) = m y K -1 (95) 

where 

m _ 1 /o°° g n exp(-2dpx - 27 / 0 T ^(z)dg)dx 
n fo° x ’ l ~ 1 exp(— 2dox — 27 Jq p(z)dz)dx 

By substituting the computed matrices ip x , rp y , R xx ( 0) and i? x „(0) into (77), we 
obtain 


^xx(w) = m y I<~' l (D-iu J M)G(iu) 
k 2(n-1) 

+ ')'rfr( II /j)GM^ T G(iw) + o(7 2 ) 

n7V j= 0 


Consequently, the spectral density matrix is given by 

#xxM = $ xx (u>) + ^(w) 

= 2dom^G(-iu})<r(r T G(iu) + 


2(n-l) 


+27 nTr^ II + 7 2 ) (96) 

j = 0 


By straightforward calculation, we know 


«r 
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Therefore, the perturbation expression of the spectral density matrix $ X xM is given 

by 


l 2(n-l) 

$xxM = G{-iu)(ra T G{iu) + 2in n [ — -( Ij)^t{G(iu)a<7 T G(iu>)} 

n7r j= o 

——G(—iu)<T<r T G{iuj)] + 0( 7 2 ) (97) 

n! 
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Chapter 6 

Infinite Dimensional Nonlinear 
Damping Models 


In this chapter, we study energy type nonlinear damping model in an infinite dimen- 
sional setting. According to the geometry of the structures considered, energy type 
nonlinear damping model is divided into two types. The following two examples are 
considered to be representatives of the two types of models, which, later on will be 
called TYPE I and TYPE H model, respectively. 

6.1 Nonlinear damping models - the formula- 
tions 

TYPE / Model: 

We consider a uniform Bernoulli beam with length L, and both ends hinged. Let 
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u(t, 5 ) denote the small deflection. Then the undamped model is given by 


ii(<,s) + a 2 uttft(t,s ) = 0 0 < s < L 

< u(i,0) = «//(<, 0) = 0 (98) 

u{t,L ) = uff(t,L) = 0 


where super-dots represent derivatives with respect to time t, and the primes deriva- 
tives with respect to s. 

We introduce the Hilbert space H = T 2 [0,X-] and the inner product defined on 
it 

[L 

[u, v] = / u(s)v(s)ds 
Jo 

Let the operator A be defined by 


Au = a 2 uW/(f, s) 


with domain 


V(A) - {u e H \ utnt € L 2 [0,Z];u(0) = u//(0) = u{L) = utt{L) = 0} 


Then, it is easy to see that 


for 


A l / 2 u = —auff(t,s) 


u € V(A 1/2 ) = {ueH \ u>r € L 2 [0,L];u(0) = u(L) = 0} 

And the total energy, the sum of strain energy and kinetic energy, is given by 

E(t) = l/2([Au,«] + ||u|| 2 ) 

= 1/2 J (a 2 u//(f,s) 2 + tl(<,s) 2 ]ds 
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If we use proportional damping A 1 ^ 2 as linear damping and E 9 (t), q > 0 as 
nonlinear damping, then the total damping is 

2^A l ^u + 7 E q {t)u 

= — 2£aurr(t,s) + {1/2 J [a 2 uft(t,sf) 2 + u^s/) 2 ]^/} 9 

Using x(t) to denote «(<,•), we can write the complete model in the following 
abstract form 

x(t) + 2M 1/a *(0 + 7 E{i)'x(t) + Ax(i) = 0 (99) 

TYPE l Model: 

We consider a flexible beam with a tip mass m at one end and with the other 
end clamped. We may assume that either the motion occurs only in the horizontal 
plane or the structure is in micro-# state. Therefore we do not have to take into 
account the effect of gravity. 

The undamped model is given by 


El 

u(f,r) H ufftf(t,r ) = 0 

P 

(100) 

El 

u///(t,£) = 0 

m 

(101) 

u(<,0) = u/(f,0) = utf(t,L) = 0 

(102) 


Notations: u(t,r) 

E 

I 

P 

m 

L 


transverse displacement at a spatial point 0 < r < L at time t > 0; 

Young’s modulus of the arm material; 

area moment of inertia of the arm material; 

density of the arm material; 

mass of the payload including the end effector; 

length of the arm material. 
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Now, in order to give a complete formulation of the problem, let us introduce 
the Hilbert space 

H = L 2 [0,L] x 1R 1 

and define the operator A by 


Aw = 

with 

V(A) ={w = 
where 


a 2 u////(-) 

—b 2 u///(L) 


for V w = 


«(■) 

U ( L ) 


ev{A) 


u(-) 

V u ( L ) 


€ H\unn(.) <e I 2 [0,T],u(0) = u/(0) = utt(L ) = 0} 


a 2 = El /p- 6 2 = El/ 


m 


The inner product on 7{ is defined by 

1 f L 


1 [L, J 

[u>i,u> 2 ] = -J / «i(r)u 2 (r)efr + — u x {L)u 2 {L) 

CL JO 0 


for 


W 3 = 


«i(0 

Uj(L) 


eH, j = 1,2 


Then, by integration by parts, it is very easy to verify that 

t L 

[Aw,w}= / |utf(r)| 2 dr > 0 Vw€V(A) 

J 0 


and [Aw,w] = 0 if and only if w = 0, i.e., the operator A is self-adjoint, positive 
definite and A~ l exists and is compact. Then by the spectral theorem of posi- 
tive self-adjoint operators with compact resolvent, there is a sequence of increasing 
eigenvalues 

u\ <u\ < • ■• 
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associated with the corresponding eigenfunctions {t/»„} such that 


Mn = 


and {V’n} form an orthonormal basis in 7f. 

With the above preparation, (100)-(102) can be rewritten as 

w(t) + Aw(t) = 0 (103) 


For linear damping part, we choose the following asymptotically proportional 
damping 

Dw = 


— auff(-) 
b 2 /auf(L) 


with domain 


V{D) = {w = 


«(•) 

u(L) 


€ H \uff € L 2 [0,l/];u(0) = 0} 


And it is easy to verify that 

r L 

[Dw, u>] = 1 /a |u/(s)| 2 ds > 0 

Jo 


i.e., D is positive definite on T>(D). 

For nonlinear damping part, we consider the energy possessed by the beam alone 

E„m = E(t)-±u(t,LY 

= 1/2 (jf |u«(l,s)|| 2 <ii + 1/a 2 jf |u((,»)| ! <li) 

Eb(t) is the energy of the beam itself, the sum of strain energy and kinetic energy. 

Then we use 7 i?j(t)iy(<) as the energy type damping. We notice that 
l/(26 2 )u(f,L) 2 is essentially the kinetic energy of the tip mass. The reason of ex- 
cluding the kinetic energy of the tip mass in damping is that internal damping 
mechanism should be dependent upon the beam material, not on any attatchments. 
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\ 


energy of the tip mass in damping is that internal damping mechanism should be 
dependent upon the beam material, not on any attatchments. 

Then the abstract form of the damped model writes 


u>(f) + 2 (Dw(t) + + Ax(t) = 0 


(104) 


Next, let us consider the basic modes and frequencies of the model (103). We 
will explicitly solve the eigenvalue problem 


Let 


A4> n = uli ) „ ip n e V(A) 


6 n (.) 

0« = | W I € V(A) 

ML) 


(105) 


then (105) is equivalent to the following two point boundary value problem (TPBVP) 

a? M f,, (r) = u>lM r ) 


b 2 <j> n r/t{L) + L> 2 n ML) = 0 

{ MO) = M{ o) = M'(L) = o 


(106) 


Let uj n = a(^) 2 . Then (106) can be rewritten as 

M"(r) = (§)Vn(r) 

+ (y )V»(i) = 0 
^n(0) = M0) = M(L) = 0 


(107) 

(108) 
(109) 


The general solution of (107) is given by 


Pn 


Pn 


Pn 


P « 


M r ) = c i,n cos -j-r + c 2(tl cosh ^-—r + c 3>n sin + c 4i „ sinh 
where c,- >n , i = 1,2, 3,4 and fl n are to be determined by (108) and (109). 
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(109) implies that 


<M r ) = ^-[(cosh ^j-r - cos —r ) - 7 „(sinh y-r - sin !j-r)\ 

where c„ is a normalizing constant and 

_ cosh 0 n + cos Pn > Q 
sinh 0 n + sin P n 

In order for <f> n to satisfy (108), /?„ must satisfy the following transcendental 
equation 


sinh P n - sin fi n - 7 „(cosh /?„ + cos p n ) 

+— ^[cosh p n — cos fi n - 7 „(sinh^„ - sin/?„)j = 0 (110) 

P L 


which is equivalent to 


TTX 3 

1 + cos fi n cosh Pn + P (cos Pn sinh /?„ — sin P n cosh /?„) = 0 

P L 


( 111 ) 


It is easy to realize that 


p n — ► tit + 7r/4 as n — ► oo 


i.e., 


<4 = A^)' — Aj )*(•» + 1/4)* ** n - CO 


Finally, the normalizing constants c n should be chosen such that HV’nll 
tedious but straightforword calculation, we obtain 


= 1. By 


c n 


= o| L + £(^)> ( i±^c^, )rj n = 12 

rn P n ' sin P n + sinh p n 


Therefore, the eigenvalue problem (105) hats been solved. 
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6.2 Some basic results 


Notations: 

E(t) = l/2([A'' 2 x(t),A l ' 2 x(t)] + [x(t),x(t)]Y, | ‘ | £He = V{A 1 ' 2 )®H 

X 

E D (t ) = l/2([Ax(t),Dx(t)] + [Dx(t),x(t)}Y I ~ ) GV(A) = V(A)<8>V(D) 


( 


T(w) = 


0 1 I u 

w = 

l 

V 


\ -7 E 9 {t)v 
Then the nonlinear damping model 


x(t) + 2 £Dx(t) + 'yE 9 (t)x(t) -f Ax(t) = 0 


can be recast into the form 


dw(t) 

dt 


= Aw(t) + r(u>(t)) 


( 112 ) 


The other notations are just as usual. 


6.2.1 Existence and uniqueness 

A. Lunardi [53] established the existence and uniqueness of solutions of the nonlinear 
infinite dimensional system 

w(t) = Cw{t) + G(iu(t)) 

by assuming 

1. C generates an analytic semigroup, which is asymptotically stable, i.e., 

sup{3?(A); A € cr(jC)} = —Uo < 0 
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2. G(0) = 0, and G/( 0) = 0; 

3. If «;(■) € C([*o,ii]; V(£)), then G(u>(-)) € C([<o,*i]; D L (9, oo)), where 
Dl(6, oo) is the interpolation space defined by 

D l {0, oo) = {x|x € X,[x] e = sup ||i , “ # £e‘ £ *|| < oo} 

0<t<l 

II x IIdx.(0,oo) = ll x ll "h [ x ]* 


The proof is essentially the apphcation of the fixed point theorem of contraction 
mapping and interpolation space theory [18] [64]. 

Therefore, to establish the existence and uniqueness of the solution of (112), it 
suffices to verify the third assumption. 

In fact, for any 


w(t) =| | € G([f 0 jti]; T)(A)) 

v{t) 


[I>(.))]* = sup ||< 1 -Mr(t)r(u;(-))|| 

0<t<l 


0 


< sup \\T(t)A I III 

o«<i \ — 7 || U j(-)|| 2, v(-) 


= 7lM•)ll 2, sup ||T(<) 
0<«<1 

< CO 


f -*> ) 

i 2(IM) ) 


Hence, 


lirw-))llD t( «.oo) = ||r(ti>(-))ii + [r(w())]» 

< OO 
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Theorem 9 For V w 0 such that ||u> 0 || < r for some r > 0, there exists a unique 
solution u>(-) € C^QO, oo); He) H CQO, oo); 'D(A)) which satisfies (112). 

6.2.2 Asymptotic Stability 

Theorem 10 Assume that T(t), the semigroup generated by A, is asymptotically 
stable, i.e. 

Iinoil < c- wt ; u! € (O,u> 0 ) 

for some u>q > 0. Then, 'iq € (0,u>o), there exists r = r(rj) > 0 such that for any 
w 0 € T>{A) with IKII < r, we have 

IKON < IKIle-"' 

Proof: For any rj € (0, w 0 ), choose ft > 0 such that rji = rj + fi < u> 0 . 

Let r = (^) 5 «+ 7 . Then, for Vw 0 € T>(A) with ||ti>o|| < r, we have 

imoii = \\T(t)w 0+1 fnt-snwis^dsw 

Jo 

< e -mt ||u>o|| + 7 /* e -ni ^ - *)||u7(s)|| 2,+2 <fs 

Jo 

Then, we have 

IH0ll e ’ ,,t < \\™o\\ + ~f f e m ’\\w(s)\\ w ds 

Jo 

< lko|| + 7lKH 29+1 / e^'IMOIM 

Then, by Gron wall’s inequality, we obtain 

IMOII ^ IKIlexpf-^i -7lKI| 2,+1 )<] 

< \\woUe-* 


□ 
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6.3 An infinite dimensional Krylov-Bogoliubov 


appr oximat ion 

6.3.1 Preliminary results 

Lemma 7 Suppose A generates an analytic semigroup T(t), which is also 
totically stable, i.e., there exists T) > 0 such that 

IPWII < e-”'; t > 0 

Then we have the following estimates for (.Dx(<), £(<)]: 

1. For linear damping model 

x(t) + 2 £Dx(t) + Ax(t) = 0 

we have 

2 £[X)i( 0 ,i( 0 ] < 1 1 «"o[ I ( 2 £(( ) J 171 , t > 0 

2. For nonlinear damping model 

x(t) + 2 £Dx(t) + 7 E q (t)x(t) + Ax(t ) = 0 

if 

A 1/2 D = DA 1/2 on V{A) 

A l/2 < pD on V(D) 

for some p > 0 , then we have 

2t[Dx(t),x(t)]< M 2 (2E(t)) 1/2 , i> 0 

where M 2 depends on £,~f,T),p and the initial state x( 0 ), x( 0 ). 


asymp- 


(113) 

(114) 
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Proof: ( 1 ) Since T(-) is analytic, we have [15] [58] 


Then, 


M, 

ll<4T(OII < -ye-*, t > 0 


A + A* 

2([Dx{t),x{t)} = -[ — - — u>(t),u>(<)] 


= — [,4u;(t), tu(t)] 

< IMHOII IKON 

= ||>tT(t)u,o|| \\w(t)\\ 

(2) To prove the nonlinear version, we consider the quantity Ei)(t) defined by 
Ed{ t) = l/ 2 ([Ax(t),I>x(t)] + [Dx(*),x(f)]) 

By (113), we can show that [ Ax,Dx] > 0 for x € P(A), and therefore 


And, we further have 

dE D {t) 

dt 


E D (t) >0, t > 0 


= -2e||Z)x(0H 2 -7^(0[^(0,i(0] 

< 0 , t > 0 


Therefore, 

E D (0)-E D (t ) = 2 £ ||T>x(s)|| 2 ds + 7 f* E 9 (s)[Dx(s),x(s)]d& 

Jo Jo 

< e d ( o) 


In particular, 


2 £ f°° ||Dx(t)|| 2 <fc < ^( 0 ) 
Jo 
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•e., Dx( ) € £ : [0,oo; H\. 

Io addition, by (113) and (114), we can have that for any v € 


iM'Ni 2 


n=0 

< p 2 52{dk,M 2 IMJI 1 
= f> 2 lii>»ll 2 


Next, we estimate 


< 

< 

< 

< 



AT(t — s)r(u»(s))<is|| 
T(t — 5)e4r(if(s))<is|| 


T(t - s^E'is) 


' *(s) 

K -2 tDx(s) 


rfs|| 


7^ 9 (0) f e^ t_J ^[||y4 1 ^ 2 i(s)|| 2 + 4 ^||£>x(5)|| 2 ] 1 / 2 ^ 
Jo 

~/E’( 0)(4{ 2 + p 2 )' 11 j‘ e-"<‘->||Di(«)||<fc 
tE’( 0)(4f* + p 2 ) 1/2 | J‘ 


+ 7^) ,/2 ^(0)4 /2 (0) 

2 77 £17 


7C0 


Therefore, we obtain 

2£[Dx(f),x(f)] = — [.4u>(/),t/;(f)] 

< \\Aw(t)\\ ||u;(t)|| 

< DWOwbll + II f* *T(t - 5)r(u>(s))rfs||]||w(t)|| 

Jo 

< [^■||w 0 ||e-" + 7 Co](2£(i)) 1 ' 2 

< M 2 {2E(t )) 1/ 2 
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which completes the proof. O 

S. Chen and R. Triggiani [31] have given the sufficient conditions on D which 
guarantee the analyticity of T(t ) generated by A. That is, if for 1/2 < a < 1, it 
holds 


PiA a <D< p 2 A a 


for some constants 0 < p\ < P 2 < oo, then T(t ) is analytic. And it is further proved 
that for 0 < q < 1/2, the semigroup is not analytic. 


Theorem 11 For the following nonlinear damping model, 

x(t) + 2 £A 1/2 x(t) + ' 1 E q {t)x{t) + Ax(t) = 0 


its solution satisfies 

x(t) = ex P (“1/ E q (s)ds)S(t)y(t) 

where S(t ) = e~ ( -'^ i , and y(t) satisfies the following undamped equation with expo- 
nentially decaying parametric excitation 

y(0 + [(1 - C)A - E q (t)A 1/ 2 + 7 2 ^(t)]y(t) = 0 (115) 

where 6(t) is a function of t which goes to zero exponentially as t —* oo, if q > 1/2. 

Proof: If we let 

x(t) = ex P(“^ E q {s)ds)S(t)y(t) 
then it is not hard to verify that y(t) satisfies (115) with 

0(t) = -l/4£ 2l (<) +?£ ! ’~ , (l)l|i(‘)ll I + %[Di(i), *(()]£>-'(<) 

Then using Lemma 1 , we have 

|0(()| < (1/4 + 2 ,)£?•(') + v'S 9 £ , - ,/! (l)[^|K||e-'" + C 0 ] 

7 1 

Therefore, if q > 1/2, |0(f)| exponentially decays to zero as long as E(t) does. □ 
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6.3.2 The Krylov-Bogoliubov approximation 

Since we are interested in the effect of nonlinear damping, we now neglect the 
linear damping part, i.e., let £ = 0. In this case, we desire to establish a Krylov- 
Bogoliubov type approximation in our infinite dimensional setting. As we know, 
Krylov-Bogoliubov approximation technique has been widely used and its accuracy 
is often satisfactory. However, it has not been generalized to multi-DOF models. 
Krylov-Bogoliubov approximation is an averaging method. In multi-DOF case, aver- 
age method does not seem to make sense because more than one natural frequencies 
exist. What we are going to do is to make use of the special form of nonlinear 
damping - energy type damping. 

We first consider TYPE I model, 

x(t) + jE g (t)i(t) + Ax(t) = 0 (116) 

By Theorem 3 , we already know that 

*(0 = a{t)y(t) 

where 

a(f) = exp(— 7/2 f E g (s)ds) 

J 0 

and y(t) solves 

y(t) + (A -I- 7 2 0(f))j/(f) = 0 

< y(0) = x(0) (117) 

y(0) = x(0) + 7 /2£*(0)*(0) 
in which $(■) is uniformly bounded, 

|0(f)| < (1/4 + 2 ? )£ J ’(i) < (1/4 + 2q)E 2q (0) 
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Next, we need to find appropriate a 0 (t) and y 0 (<) such that x(t) is approximated 


by 


£o(0 = ao{t)y 0 (t) 


First, since a{i ) is slowly varying, we make our first approximation 


i(t) « a(t)y(<) 


Consequently, 

i(0) = y(0) (118) 

Next, we make our second approximation by letting y(t) « y 0 (t), where yo(t) is 
obtained by dropping 7 2 0(t) in (117) and replacing the second initial condition in 
(117) by (118). That is, j/o satisfies 

yo(t) + Ay 0 (t) = 0 

< yo(0) = * o (0) (119) 

jfo(O) = x(0) 

Based upon the above two approximations we made, we obtain 

m = « 2 (<)/2([^o(0 5 yoW] + l|yo(0ll 2 ) 

= «(<) 2 

= “'mm 

where 



Then, by definition, 

a(t) = exp(-7/2£ , (0) f a 29 (s)ds) (120) 

JO 
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We naturally let the solution of the above integral equation be a 0 (t). 

Therefore, differentiating (120), we realize that oo(<) solves the following initial 
value problem 


«o(0 = -i£'(o)a2’ +, (() 

«o(0) = 1 


( 121 ) 


One can easily find the solution as 


<*,(<) = (1 + 


Therefore, the Krylov- Bogoliubov type approximation of TYPE I model is given 

by 

*o(0 = (1 +79^ 9 (0)t)"^yo(i) 

From this Krylov-Bogoliubov approximation, we can see, without linear damp- 
ing, the free response of TYPE 1 model still goes to zero, and the decay is in the 
order of t~*i. Furthermore, larger value of q implies slower decay. 

For TYPE E model, similar result can be obtained. 

In fact, through the same procedure as above, we cam see that the yo{t) we choose 
should also satisfy (119) in this case, while the integral equation for a 0 (t) is slightly 
different. 

Under the same approximations, we have 

EM = E(t)-±i?(t,L) 

= aS(i)[£(0)-ii5(M)] 

where 

Sto(<) = f “° (< ’ s) ) 0<s<I 

V / 
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Then, the integral equation for a 0 (<) is given by 

a 0 (t) = exp(- iftf El(v)dtt) 

= exp|-7/2j f o al'(tl)(E(0) - £))’<«/] 

from which, one can easily find, 

ao(t) = [1 + 79^ (£(0) - -^i%(tt,L)) q dti)-b 

Thus, the Krylov-Bogoliubov approximation for TYPE H model is obtained. 

Notice that since uo {t-L) is periodic in t and 

E(0)-^u 2 o (t,L) = l/2\\Y 0 m 2 E-^l^L) 

= l/2[Ay 0 (t),y 0 (t)] + J q \u 0 (t,s)\ 2 ds 

> 0 

we can realize that as t ► oo, 

jlm)--^ii(V,L)]HV = 0(t) 

Therefore, the Krylov-Bogoliubov approximation indicates that the free response of 
TYPE II model also goes to zero at the rate of 0(t ~ ^») as t — ► oo. 

In the above, we have used the simple fact that if /(<) > 0 and is periodic with 
period T, then 

6.3.3 The Error Estimate 

Next, the natural question we ask is what the error of the Krylov-Bogoliubov ap- 
proximation is. For simplicity, we only study the error in TYPE I model. 
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Before the estimation of the error, we first introduce a correction term yo in the 
initial condition of (119) so that yo(t) now satisfies 

y Q (t) + Ay 0 (t) = 0 

i y o (0) = x(0) (122) 

y o (0) = x(0) + y 0 

The purpose of introducing the correction term yo in the initial condition will be 
explained later. 

Now our result is 

Theorem 12 For TYPE I model, the error of the Krylov-Bogoliuhov approximation 
is given by 

||x(<) - XoWII < o(()II-4‘ 1/! (7/2^(0)x( 0) - Vo) 1 1 + 7 2 Bo((,l) (123) 

where 

limS 0 (f,7) < —E 2q+1/2 ( 0)[^ 4 + 29 + q/2]t 

-•y— *0 U?X 

The error contains two parts, one part is of the form C'fH. Therefore if f < I/7, 
then this part of the error is still of the order of 0(7). The other part of the error, 
a(i)\\A~'l 2 (i l2E q (Q)x(ti) — yo)||> is decaying with time because a(t) is. In addition, 
by choosing 

V„ = y/2E*(0)x(0) + 0(7 J ) 

we can keep the first part of the error in the order of 0( j 2 ), hence the total error is 
of the form 0( y 2 ) x (t + 1). 

To further understand the reason of introducing the correction term yo, recall 
that the first approximation we made in the derivation of the Krylov- Bogoliubov 
approximation is 

i(t) = a(t)y(t) 
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since a(f) is slowly varying. Based upon this approximation, the initial condition 


y(0) = x(0) + 7 /2£ ? (0)x(0) 

is replaced by 

y o (0) = x(0) 

This change of intial condition, on one hand, made the following derivation possible, 
on the other hand, introduces an error (in the order of 0(7)). 

Through error estimate, we realize that this error can be minimized by choosing 

y 0 = 7 /2£ 9 (0)*(0) 

Therefore, the correction term yo plays the role of compensating for the error intro- 
duced by the “slow varying” approximation. 

PROOF: First we estimate ||y(f) — J/o(OII- We know that 

y(t) = C(t)x(0) + 5(0(i(0) + 7 /2^(0)x(0 ))- 7 2 f e{r)S{T)y{r)dr 

Jo 

y 0 {t) = C(t)x(0) + S(f)(x(0) + y 0 ) 

where C(t ) : H — ► H and S(t ) : H — ► T>{A X ^) are cosine and sine operators 
defined by, Vx € H, 

OO 

C(t)x = J 2 COSUJ n f[x,V> n ]V>„ 

n=l 

S{t)x = Sma;rtj! [x,V>n]^n 

n=l U>n 

Obviously, we have 

||S(<)x|| < P" 1/2 x|| for xeH 
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Then, 

IM<) - soMII 

< ||S(i)( 7 / 2 E , ( 0 )x( 0 )-ji„)|| + 7 J (l/4+2 g )£; J ’(0) £ ||S(r)y(r)||<ir 

< ||A- 1 ' j ( 7 /2£:>(0)x( 0) - Mil + 7 2 (l/4 + 2q)E^0)J‘ ||A-'/ 2 y(r)||^124) 


Secondly, we estimate |a(<) — ao(f)|. For this purpose, we consider 



is a unitary operator from He to He- 

From the above equation, we can easily know that 

l|K(0lll = 2E(0) + 7 ^( 0 )[x( 0 ), x(0)] + 0 ( 7 2 ) (125) 

In addition, from 

*(0 = -~1 l‘lE q (t)a(t)y{t) + a(t)y(t) 

we know that 

l|i(<)H 2 = ° 2 (<)lli(0 - 7/2£ , (0»(0ll J 

= « 3 (i){lliKi)lP - 7^’(t)[y(<), »(«)] + 7 2 /4-E 2 »(*)lly(0ll 1 } (126) 

Then, by (125) and (126), we can rewrite the energy corresponding to TYPE 1 
model (116) as 

E(t) = l/2{[Ax(0,x(0] + ||x(*)|| 2 } 
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= i/ 2 « ! (i){[/is(<). »(«)] + llsKOII 2 - tB’WIkW. »(<)] + 7 74£ S> (0IM0II 2 } 
= a J (iMi/2||y(()lll - 7 /2£ J (0M0.y(0] + 0( 7 ! )} 

= a ! (i)[£(0)+ 7 /(() + O(7 2 )] 

where 

m = 1/2{£«(0)[X(0), x(0)] - £'(0M0, «<)]} 

Then, a(t) can be rewritten as 

a(t ) = exp(— 7/2 J E q {r)dT) 

= exp {- 7 /2 J‘ S<(t)[E(0) + 7 /(t) + 0(l 2 )]’dr} 

from which, a(t) can be easily solved to be 

a(t) = [I+7 q f [E( 0 ) + 7/(7-) + 0(7 2 )] 9 dr] - ^ 

Jo 

= [1 + ~rqE q (0)t + 'y 2 q 2 E 9 ~ 1 (0) f f(r)dr + 0( 7 3 )]'^ 

Jo 

Therefore, we obtain 

WO - «o(OI < 7 2 ?/2£ ? -'(0)[jf ' f(r)dr + 0( 7 ) ) (127) 

Finally, using (124) and (127), we obtain the error estimate as 

I WO - *o(OII < «(0IW0-»»(0ll + W0-«o(0IIM0ll 
< a(()||/r‘' 2 ( 7 / 2 £«( 0 )x( 0 )-!i„)|| 

+7 2 (l/4 + 2q)a(t)E 2 *(0) £ \\A-'l 2 y(r)\\dr 
Wq/2E<-\0)[£ |/(r)|dr + 0( 7 )] || Vo (<)l| 

= «(0||il- ,/, ( 7 /2£'(0)*(0) - jb)ll + 7 2 Bo((,7) 
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where 


B„(i, 7) = (1/4 + 2q)a(t)E 2 '(0) J* \\A~ l ^ 2 y(T)\\dT 

+q/2E'-'(0)[J‘ |/(r)Mr + 0( 7 )]|Mi)ll 

By noticing that for 7 — * 0, 

|/(i)| < l/2£>(0)(||x(0)||||i(0)|| + |M()||||»o(()||) 

< ^(IM ,/2 x(o)|| ||i(0)|| + |M , 'V1)|| lljb(OII) 

ZU3\ 

^ £’(0) ||A 1 /2x(0)|| 2 + ||x(0)|| 2 IIAV^IP + llt/oW, 

- 1 2 2 ’ 

< — E 9+1 (0) 

Wl 

we can obtain by easy calculation that 

lim B 0 (<,7) < —[9/2 + 1/4 + 2, )£ wl/2 (0)i 

*r-»0 uj 1 oil 

□ 

6.4 Frequency response - single frequency exci- 
tation case 

We consider the following infinite dimensional nonlinear damping model 

x(<) + 2f A'/’x(t) + 7(H*(0,*W1 + l|i(<)ll 2 )'i(0 + 4lx(<) = E(t) (128) 

For the reason of simplicity, we consider strictly proportional linear damping and 
total energy ( TYPE I) nonlinear damping so that we can easily obtain its modal 
decomposition. 
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We are interested in the stationary amplitude response of each mode to a single 
sinusoidal input to the system. To specify the excitation input E(t), let B : 2R 1 — ♦ 
H be a linear bounded operator. Then we let the excitation be of the form 

E(t ) = tB cos u>t 

where t > 0 is a small parameter in the order of ( and 7. That is, we let 

( = (of, 7 = 7o c 

\ We order the input amplitude as 0(e) for the following reasons: 

1. For a weakly damped system, a small amplitude (0(e)) excitation at the nat- 
ural frequency produces a relatively large (0(1)) amplitude response. 

2. In the actural system, large oscillation are limited by damping. Thus to obtain 
a uniformly valid approximate solution of this problem, we need to order 
the external excitation amplitude such that, in the perturbation procedure, 
excitation term appears in the same equation as damping terms. 

Therefore, we are considering the stationary response of 

x(t) + e[2f 0 A}! 2 x(t) + 7 0 £ , (t)i(f)] + Ax(t) = eB cosut (129) 

The method we will use for the following analysis is multiple scale method. Define 

T 0 = t\ Ti = et 

and the notations 

Aj/fTi, = ^ D.HTo.Tt) = ^ 
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Let the solution of (129) be 


x(t) = Xo(7o, Xi) + exi(To, Tj) + • • • 

then 

x(t) = D 0 x 0 + e(DiX 0 -+ D 0 xi) + 0(c 2 ) (130) 

x(t) = DqX o + e(2D 0 Dix 0 + D^x i) + 0(e 2 ) (131) 

Then substituting (130)-(131) into (129) and equating those terms withe same order, 
we obtain 

e° : DqXo(To, Xi) + Aro(Xo, Xi) = 0 

e 1 : Z)qIi(Xo, X i) -(- v4xi(Xo, Xi) 

= —2D 0 Dix 0 - 2( Q A l/2 D 0 x 0 - 7o([>lxo,xo] + IIAjXoIH’A^o + B cosuit 

Obviously, Xo(X 0 ,Xj) and x^XojXx) can be written in the form 

OO 

£o(Xo,Xi) = Y Qm(Xo, Xl)V> m 

m = l 

»i(r„r,) = 

m = l 

Then, {a„(7o,Ti); n = 1,2, — } and {6„(X 0 ,Xi); n = 1,2,--} satisfy 

X*o a n(Xo, Xj) + w 2 a n (Xo, Xj) = 0 
DlKWoW + tfMToiTx) 

OO 

= -2D 0 X>ia„ - 2(oui n D 0 a n - 7o[^ + (D 0 a m ) 2 ] q D 0 a n + cos ut 

m—\ 

n = 1,2, •• • 

Obviously, for n = 1,2,-**, 

a n (Xo,X x ) = A n (Ti)e iu>nTo + A n (T l )e~ iu ' T ' 
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from which, we know 


E^a^(r„,r 1 ) + [n„<. m (r„,T 1 )] 1 

m=l 

oo 

= £«J,|2.4»(7i)|» (132) 

m= 1 

Then the equations for 6 n (!To,Ti)’s become 

DX(To,T,)+^MTo,T,) 

= e-“" T "[-2iu, n /l n /(r 1 ) - 2,(<vlA„(T,) - 7o(E ^|2A„| 2 )’iu.„A n (r,)] 

m=l 

+B m tl> n / 2e iwt + cc (133) 

where cc stands for complex conjugate of the preceding terms. 

Next, we assume that the excitation frequency u> is close to a particular natural 
frequency, say u>j, i.e., without loss of generality. 


U) = U)\ + to 


where a is a detuning parameter, quantitatively describing the nearness of u> to u>i. 

Then bi(To,T\) satisfies 

^6 1 (T 0 ,r 1 ) + w 1 2 6 1 (r 0 ,r 1 ) 

= (134) 

m=l 

+B m ^/2e iaTi ] + cc (135) 

If the coefficient of e ,u ' 1 T ° is not zero, then the right hand side of (134) is 
of the form /(J\) cosu>i7o, which would produce terms like T^ 1 cos u>i To and/or 
Tq 1 sinwjTo. Such terms are called secular terms , which are obviously unbounded. 
Since a positively damped system cannot have unbounded response corresponding to 
a bounded input, secular terms should not appear. Therefore, in order to eliminate 
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secular terms, we need to choose A^Ti) such that the coefficient of e“ 1 ' lTo vanishes, 
i.e., 

—2iu>\A\i(T\) — 2i(ou%Ax(Ti) — 7o( ^ ^^,|2i4 m | 2 ) , zu;ii4i(Ti) 

m=l 

+B*rl>\/2e i<,Tl = 0 (136) 

For the same reason, for n > 2, j4 n (7i) should satisfy, noticing + «r is 

away from u„, 

- 2iu n A n t(Ti) - 2ii 0 ulA n {T l ) - 7o(f ] u^A^yi^A^) = 0 (137) 

m=l 

and hence, for n > 2, b n (T 0 ,Ti) satisfy 

DX(To,T,) + u*MT 0 ,T t ) = B'tp n coswt (138) 

Next, in order to solve {A n (Ti); n = 1,2, - - • } from (136) and (137), let 
A n (T,) = ^Pn(T,)e iMT ‘\ n = 1,2, 

Then substituting A n {Ti) into (136) and (137), collecting the real and the imaginary 
parts, one obtain the following systems of differential equations for 

{/>»(ri), A(7i); « = 1,2,...} 

^ = -{ou>,p,(r I ) - 7 o/ 2(£» =1 ^^(roiVifr.) 

+‘^-sin-)’i(r,) (139) 

MT))*^ = <TPi(ri) + ^cos7,(r,) 

= -fc^nOT.) - 7o/2(E“ = ,^A(r.))Vn(ri) 

(140) 

/>n(r,)i^ = erp n (7j) n = 2,3,-- 
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It is obvious that for each solution of (140), it holds 


lim p n {T\) = 0, n > 2 

T\ — *oo 


For (139), by setting ^ = 0 and = 0, we know that at steady state, pi 

and 7i satisfy 


0 = -£o^i?i - 7 o/ 2 (u>! pif q pi + cos 7 a 

0 = ap, + Blix cos 7, 


(141) 


Then, eliminating 71 from (141), we obtain the equation for the steady state har- 
monic amplitude of the first mode 


pH* 2 + [{o“>i + 7o/2(wi/>i) 2, ] 2 J = (^) 2 (142) 

This is called frequency response equation. It is easy to see that the frequency 
response equation has a unique positive solution p\ for each value of cr, denoted by 
g{o\uj\, The plot of p\ in terms of o is called frequency response curve. From 

(142) we can see that p\(cr) = <7(er|u>i, B m x^i) is an even function of a. 

In fact, the existence and uniqueness of a positive solution of (142) becomes 
trivial if one rewrites (142) in the form 

2 _ [B*1> i/(2o>i)] 2 

1 0-2 + [£okh + 7o/2(wa/)x) 2 «] 2 

Next, let us examine the steady state response of the other modes. We already 

know that p n = 0 for n > 2 at steady state. But what about the 0(e) order term? 

Solving 6„(T 0 ,Ti) from (138) gives 



b n (T 0 , T\) = C n (7i) co S (u> n 7o + MTa)) + 


u > 2 — uS 1 


COS uTo 
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oo. 


Then it is not hard to realize that C n (T\) — ► 0 as 7\ — ► 

Therefore, for u> = u>i + ecr, the nth (n > 2) mode stationary amplitude is of 
order 0(e), given by 


Pn = e 


|£>«| 




u> 2 | 


+ 0(e 2 ) 


Similarly, we can infer that if the excitation frequency u> is away from u>\, then 
the steady state amplitude of the first mode is given by 

*(*■>)-< ■ tPQ ,. +0(t 2 ) 

M-u> 2 \ 


In summary, the frequency response curve of any mode, corresponding to a single 
frequency ui excitation is given by 

r 

g((j\u > n , B‘ij) n ) when u; = u> n + to 
p n (u) = < (144) 

when |a>-u? n | = 0(1) 


6.5 Frequency response - multi-frequency exci 
tation case 


Next, we consider the case in which the external excitation E(i) contains M distinct 
frequencies, and each of them is close to a natural frequency. Say, they are u n + 
ecr n , n = 1, 2, • • • M. Then the excitation takes the form 

M 

E(t) = e f n Bcos(u n t + ea n t + r n ) 

n= 1 

where /„ and r n are certain constants. 

Through similar procedures and with the same notations as in the single fre- 
quency excitation case, we can obtain 

Dlb n (T 0 ,Ti) + u%b n (T 0 ,Ti) 
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C-3L 



_ e«"" r °[_2iu> n ,4 n /(Ti) - 2i( 0 U>lA n (T 1 ) - 7o(S ^ 2 |2A m | 2 ) 9 tU; n A„] 


m=l 


M 


+ B' Vn/2 53 / J V w > To+ *' r,+, ' T > n = 1, 2, • • • 

j=i 

In order to eliminate secular terms, A n (Tj) should satisfy 
— 2lU> nJ 4 n /(7i) - 2i^ 0 ^i4 n (T 1 ) - 7o(5Z 

m=l 

+£> n /2/ n e ,(£, " Tl+T "> = 0 

for n = 1,2, • • • , M 

-2zw„i4„/(7\) - 2^ 0 w*A n (r 1 ) - 7o(53 w 2 |2A m | 2 ) ? tu> n A n (:Ti) = 0 

m=l 

for n = A/ + l,M + 2,--- 


Similarly, by letting 


A n (T 1 ) = p n (T 1 )/ 2e*'^ T >) n = 1,2, ■ 


we can know that 


and for 1 < n < M, 


T\ 


lim p n {T\) = 0 for n > M + 1 


= -W,(r,)-r./2(E^A(r,))v,(T,) 

* m= 1 

+B m tp n / (2 u>„)/ n sin 7 „(Ti ) 

^(r,)^p = ff.A.(r 1 ) + B-A,/(2u W )/.a»7„(r 1 ) 

where 7„(7 1 1 ) = <r„ri + r n - 

Again, after reaching stationarity, />„, 1 < n < M satisfy 

<({<«, + W2Q>M>*l 2 + <*>VI = 

j=l 


(145) 


88 



Therefore, in the case of multi-frequency excitation, the stationary harmonic 
response amplitudes p n , 1 < n < M satisfy the above coupled nonlinear equation. 
In order to study the solutions of (145), we introduce the notation 

G = 'ro/2(Y^^jPj) v 

j = i 

Then, from (145), we obviously have 


2 2 (fnB^n/ 2? 

“nPn = 


n = 1,2, • • • ,M 


( 146 ) 


(£o^n + G ) 2 + 

On both sides of (146). summing up from 1 to M, raising to the power q and 
multiplying by 7o/2, one obtains 


G = 70/2E J ' 


(147) 


■a (6 m. +g) j +^' 

It is easy to see that for each a 2 = (a 2 , o\, • • • , <Tm), (147) has a unique positive 
solution denoted by G(a 2 ). Therefore, for 1 < n < M, p n satisfies 

n /( 2 u , n )) 2 


Pl(n = 


(148) 


ff2 + Kow B + G(5»)] a 
Pn(° 2 ) is single valued because of the uniqueness of G^a 2 ) for each a 2 . 

To study the behavior of p n (B 2 ), we first notice that by differentiating (147) with 
respect to a 2 , we realize that 

dG(a 2 ) 


da 2 


< 0 j = 1,2 , -,M 


i.e., G(B 2 ) is decreasing as a 2 increases with the other a*’s fixed. 

Next, we restrict our attention to the case M = 2, the interation of the first two 
modes. As a 2 increases, G(ai , a 2 ) decreases, and hence Pi{o\,o 2 ) increases. As a\ 
increases, p 2 (a j,cr 2 ) increases because G{<7\,o 2 ) decreases. Then it must be the case 
that />i((7i,<t 2 ) decreases, since G{&\, a 2 ) decreases. Therefore the plot of p\(<T\,a 2 ) 
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is a saddle surface. For fixed <7 2 , the curve Pi(*, 02 ) is still bell shaped, reaching 
maximum at <7\ = 0. As a\ decreases, the energy possessed by the second mode 
becomes bigger which results in bigger damping in the nonlinear energy type damp- 
ing, therefore due to coupling the first mode amplitude becomes smaller. Similarly, 
p 2 (<7i, 02 ) i s a l so a saddle surface. 

Numerical Results 

For any fixed n, 1 < n < M, in order to find the frequency response p n {B 2 ) 
from (148), one has to resort to numerical methods to find G(B 2 ) first for each B 2 . 
Here, we pointed out that if we use fixed point iteration starting from 0, then we 
have monotone convergence, hence avoiding the common phenomenon in solving 
nonlinear algebraic equation that the convergence depends on the choice of initial 
data. 

Specifically, we define 


M 


r(<r 2 ;x) = 7 o/2[J3 


n= 1 


+ (£oU>„ + x)2 J 


and 

| r 0 (<7 2 ) = r(<7 2 ; 0 ) 

I T n (B 2 ) = r(or 2 ; r n _i(^ 2 )), n = 1,2, • • • 

It can be shown that r 2n (<? 2 ) is monotone decreasing as n increases, and 
r 2 n+i (^ 2 ) is monotone increasing as n increases. In addition, it holds 


r 2n+ i (<? 2 ) < G(B 2 ) < T 2n (B 2 ), n = 0, 1 , 2, • . • (149) 


In fact, first it is obvious that 


r,(«7 2 ) < g(b 2 ) < T 0 (<? 2 ) 
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Then, we can easily show that 


G(s>) < r J (5 2 ) = r(^r J (^))<r(s 2 ;0) = r„(5 J ) 
G(<? 2 ) > r3(5 2 ) = r(5 2 ;r 2 (5 2 ))>r(5 2 ;r 0 (? 2 )) = r,(5 2 ) 

Then, by assuming 

G(5 2 ) < r Jn (5 2 ) < r 2n .,(5 2 ) 
g(j 2 ) > r, n+1 (s 2 ) > r 2 „_,(o 2 ) 


we obtain 

G(? 2 ) < r 2n+2 (5 2 ) = r(j 2 ; r 2n+I (5 2 ))<r(s 2 ; r 2 „_ 1 (5 2 )) = r 2n (? 2 ) 

G(? 2 ) > r !n+3 (s 2 ;r 2 „ +2 (? 2 )) > r(ff 2 ;r 2 „(5 2 )) = r 2 „+,(? 2 ) 

Therefore, by induction, we have shown that {T^o 2 )} ({f' 2 n +i{o 2 )}) is mono- 
tone decreasing (increasing), and (149) holds. 

A few numerical examples are made, in which we have chosen M = 3, q = 
3/2, 70 = 2, f n B m ipn = 2, for n = 1,2,3. From the computer tests we realize that 
when £o is not small, the convergence is quite satisfactory, while the convergence is 
very slow for small £o- Therefore, in choosing e, we should let e = £/k where k > 2 
to guarantee fast convergence. 
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CM 

II 

O 

6 = 2 

6 = 1 

It 

0 


0\ — 0 

a x = 1 

<7l = 1 

0\ — 3 


O 

II 

CM 

b 

ct 2 = 1 

<r 2 = 2 

CM 

II 

CM 

b 


0 

II 

CO 

b 

<7 3 = 1 

<73 = 3 

<7 3 = 1 

Solution 

0.333 

0.236 

0.448 

0.848 

To 

0.662 

0.323 

0.811 

1.588 

r. 

0.195 

0.211 

0.278 

0.379 

r 2 

0.433 

0.244 

0.563 

1.356 

r 3 

0.280 

0.234 

0.384 

0.483 

r 4 

0.368 

0.237 

0.488 

1.241 

r 5 

0.314 

0.236 

0.424 

0.547 

r 6 

0.346 

0.236 

0.462 

1.169 

r 7 

0.326 

0.236 

0.439 

0.592 

r 8 

0.338 

0.236 

0.453 

1.118 

r 9 

0.331 

0.236 

0.445 

0.627 


By ploting p\(cr\ , (^l), we verified that it is indeed a saddle surface. However, 
the increase in er 2 direction for fixed cr : is very small, which is also clear from the 
corresponding contour plot. The decrease of Pi(&i , direction is much more 

significant than that in <72 direction. This indicates that even though there exists 
coupling due to nonlinear damping, the coupling effect is rather weak. 
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Chapter 7 

Active Damping of Flexible 


Structures via Saturating 
Acturators 

7.1 Introduction 

In the previous chapters we have studied nonlinear passive damping problem. In 
applications, active dampers are often used to enhance the stability of a flexible 
structures. Moreover, due to actuator saturations, active damping become non- 
linear. Therefore, it is also important to study nonlinear active damping, as well 
as passive damping. In this chapter, we study the active damping aspect of the 
SCOLE problem [4] - a recent NASA project. The primary objective of the SCOLE 
(Spacecraft Control Laboratory Experiment) problem includes the task of directing 
the line-of-sight of the shuttle/antenna configuration towards a fixed target (Figure 
1). Due to the facts of very small passive damping in the supporting truss structure 
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and the micro-5 environment in the orbit, structure vibration is inevitable after each 
slewing maneuver. In order to maintain the prescribed pointing accuracry of the 
antenna line-of-sight, active dampers are required to enhance the structure stabil- 
ity. Only active damping problem will be studied and colocated sensor/actuator 
arrangement will be used. The controls are force and moments actuators and the 
sensors are rate gyros. 

A distributed parameter model is used in this investigation. In Section 2, the 
continuum model and problem formulation are presented. 

In Section 3, a group of rather weak sufficient conditions for strong stabilizability 
is presented. Although other sufficient conditions for strong stabilizability has been 
obtained in [7], the present conditions are much weaker in the sense that internal 
damping is no longer required to be positive definite, and, in fact, can even be zero 
while strong stabilizability can still be achieved by active damping. 

In order to understand the nature of active damping, in particular, the nature of 
saturation type active damping, the mode excitation problem is studied in Section 
4. In this study, some notions in classical feedback control such as Characteristic 
Equation and Root Locus axe extended to our distributed parameter system, which 
is a feature of this work. As we will see, the root locus can provide an insight of 
the nature of active damping. The root locus method has been a powerful and 
useful approach for the analysis and design of finite dimensional control systems. 
However, the notions of Characteristic equation and its Root locus have not been 
extended to the study of distributed parameter systems, due to the difficulty caused 
by its infinite dimensional nature. In this work, we have taken the advantage of the 
fact that, although the system is infinite dimensional, the number of actuators and 
sensors used is always finite (here, we have excluded those applications in which 
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Figure 1 . Shuttle Orbiter/ Antenna Configuration 


distributed actuators and/or distributed sensors are used). 

While much work has been done in the active damping of distributed parameter 
oscillation systems by using linear feedback control [18], [39], [53], from practical point 
of view, actuators can be linear only in certain range (small amplitude) and become 
saturated for large amplitude. Therefore, we need to take the nonlinear (saturating) 
nature of the actuators into consideration in control design. Generally, actuator 
saturation brings many difficulties. Nowadays, such systems are still designed by 
intuition, experience, and simulation using trial and error. Their effects on the loop 
response are still poorly undrstood from a theoretical point of view. In Section 4, 
the effect of actuator saturation is studied and compared with linear damping case. 

In Section 5, the notions of Characteristic Equation and its Root Locus are ex- 
tended to general multi-actuators/sensors case. The generalization is based upon 
an operator inverse identity. 

Since sensor noise is inevitable for most sensors, in Section 6, we consider sensor 
noise problem in active damping through direct connection. By studying the effect 
of sensor noise on the steady state antenna motion, a design guideline is provided 
for the choice of the feedback gain constant in active damping. 

Section 7 summarizes the conclusions of this work. 

7.2 A Simple Continuum Model 

The model we will consider is based upon a simplified version of the SCOLE problem 
- a flexible truss, clamped at one end, with an offset antenna at the other end, has 
only bending motion in a plane. The truss structure is modelled by an equivalent 
uniform Bernoulli beam of length L along z-axis, extending from r = 0, the clamped 
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end, to r = Z, the antenna end. The antenna is considered as an attached tip mass 
with mass m. With r, 0 < r < L, denoting the spatial variable along the 2 -axis 
and t denoting time, let u(t,r) denote the displacement of the truss in y — z plane. 
For simplicity, we suppose we use only one control force actuator which is located 
at the antenna end. Then the beam deflection u(f,r ) sloves the following boundary 
coupled linear partial differential equations: 

FT 

u(i,r) + —u""(t,r) = 0 

p 

FT 

u(t,L)-—u'”(t,L) + z(t) = 0 
m 

tt(t,0) = 0 (150) 


t/(t,0) = 0 


«"(*,£) = 0 


where super-dots represent derivatives with respect to time t, and the primes deriva- 
tives with respect to r. Here we use z(t) to denote control force and E, 7, p are 
Young’s modulus, moment of inertia and density of the beam material, respectively. 

In order to study this problem systematically and rigorously, we need to refor- 
mulate this problem in a Hilbert space setting. 

First we introduce a Hilbert space 

H = Z 2 [0,Z] 0 2R 1 


and the inner product on 77, 

1 rL 1 

[x l5 x 2 ] = J o (r)u 2 (r)dr + — C!C 2 


where 



€ 77, j= 1,2 


97 



and 


a 2 = El /p, b 2 = Elf 


m 


For a general element x 6 H, its scalar component does not have to have any relation 
with the Z 2 [0, Z]-function component. 

And, we define the stiffness operator A by 


Ax = 


aV"'() 

-b 2 u"'(L) 


for 


with 


V(A) = {x = 


( 


u(-) 


\ U (L) 

Next, we define the control operator B as 


u(-) 

Vi = | I <E V{A) 

u(L) 


€ H I e L 2 [0,L], «(0) = u'(0) = u"(L) = 0 } 


B . IB} — ► H, Bz = 


, zeIR 1 


i.e., the control is applied only to the beam tip. Obviously, Af(B) = {0}. 
Under the above notations, (150) can be written as 


i(f) + Bz(t) + Ax(t) = 0 (151) 

It can be shown that A is self-adjoint, positive definite on T>(A) and A -1 exists 
and is compact. Then by the spectral theorem of positive self-adjoint operators with 
compact resolvent, there is a sequence of eigenvalues of A (natural frequencies) 

0 < w 2 < < • • • — ► oo 
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associated with the corresponding eigenvectors {V\», n = 1,2, •••} such that 

Aip n = ulxj> n , n = 1,2, •• • (152) 

and, furthermore, {ip n , n = 1, 2, • • •} form an orthonormal basis in H. 

The eigenvalue problem (152) can be solved to give 


in which 


V’n 


l Mr) ) 

( ML) ) 

a(^) 2 , n = 1,2, •• • 


, , , 1 r t. @ n @ n ( • k @ \1 

<(>n( r ) = — [cosh — r — cos — r — 7 n (smh — r — sin — r)J 

Cyx Li Lj Li L/ 

and {/?„, n = 1,2, • • •} are the solutions of 

ttiQ 

1 + cos fl n cosh /?„ H r— (cos /?„ sinh /?„ — sin /?„ cosh /?„) = 0 

pL 


In fact, we also have 


\P n — (n — 2 4- l/4)?r| — ► 0, asn-too 


In the above, the constants c„, 7 „ are defined by 

«w = 4 L + U±-?XIY» 

m p n 

_ cosh fin + COS 0 n 
7n ~~ sinh P n + sin p n 
x _ 1 + cos ftn cosh 

- sinh Pn + sin p n 

Later on, we will need the following realtions which are not difficult to verify, 


ML) 


2 pLXn 1 

men Pn 
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= 0 ( tJ 


= <4) ! 


as n — ► oo 


p! m % 


(153) 


b 2 


a"' 1 /L + p/m$ 


• 0( k ] 


as n — ► oo 


(154) 


Sometimes, we need to write (151) in the form of first order system 


' *(0 \ _ 

K *(<) J 

The underlying Hilbert space is 


d_ 

dt 


0 / 
—A 0 


x(f) 

x(f) 


+ 


-B 


z(t) 


He = T>(A 1/2 ) ® H 


equipped with the inner product 


[u>i,u> 2 ]e = [j4 1/2 xi,j4 1/2 x 2 ] + [yi,y 2 ], v>j = I 1 | € H E , j 

\ yj 

For later convenience, we define the operators 


= 1,2 


/ 


A = 


0 I 
-A 0 


, with V(A) = V(A) <g> H 


and 


B = | ° \ : IR 1 —* He 

-B 


7.3 Strong Stabilizability Conditions 

In this section, we present a group of sufficient conditions for the strong stabilizabil- 
ity of general distributed parameter oscillation systems, not only for the particular 
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model presented in the last section. We concentrate on the following abstract wave 
equation version of an infinite dimensional linear oscillation system 

x(t) + Dx(t) + Ax(t) = Bz(t ) z eJR m (155) 

on a Hilber space H. Last section can be considered as a simple example of reducing 
a concrete PDE model to the above abstract wave equation model. 

In (155), A is a stiffness operator with domain V{A), and is generally nonnega- 
tive definite with compact resolvent. There exist a sequence of natural frequencies 
{u;„, n = 1,2, •••} and the corresponding linear natural modes {rp n , n — 1>2, •••} 
such that 

1. Axfrn = u*xl> n n = 1,2,- • 

2. < u >2 < • • •, and lim n _ 00 w„ = oo; 

3. {V’m n = 1,2, • • •} form an orthonormal basis on H. 

From now on, we assume u>i > 0, i.e., there is no rigid body mode. 

D is the linear internal damping operator, which is nonnegative definite on its 
domain V(D). B : 2R” 1 — ► H is a finite dimensional linear operator. 

The energy of the system is defined by 

£(()=l/2[||A , ' J I (<)ll I +l|i(<)ll J l 

The question we want to answer is: taking the actuator saturation into consid- 
eration, what kind of feedback stabilizing control we should use, in order to make 
the closed-loop system strongly stable. By strong stability, we mean, for any given 
initial data, (x(0), i(0)), the corresponding closed-loop system response satisfy 

lim £(t) = 0 

<— oo v 
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Our main result of this section is the following 
Theorem 13 Under the following assumptions: 

1. /(•) is Lipschitz continuous; 

2. [/( x )» x ]r"’ > 0, x ^ 0; 

3. M{ {D + BB-) |e„) = {0} 1 

the following rate feedback with colocated sensors/actuators 

z(t) = ()) 

strongly stabilizes (155), i.e., the closed-loop system 

x(t) + Dx{t ) + Bf(B m x(t)) + Ax{t) = 0 (156) 

is strongly stable. 

Here E n denotes the subspace spanned by those natural modes %l> nj corresponding 
to the natural frequency u>„ . 

In particular, if all ui n ’s are distinct, then Assumption 8 can be replaced by 
Assumption 3' : [DxJ> n , ip n ] + ll-fiV'nll 2 >0 n = 1,2, • ■ • 

Before the proof, we make the following remarks: 

1. In this theorem, we do not require the passive (internal) damping D to be 
positive definite. In fact, even we neglect the internal damping ( D = 0) in the 
modeling, (156) is still strongly stable as long as 

||B*^ n ||>0, 1,2,.-. 

and Assumptions 1, 2 are satisfied. 

l M(P) stands for the null space of the operator P, and (D + BB’)\e. stands for the restriction 
of the operator D + BB * on the subspace E n . 
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2. Assumption 3 is also a necessary condition. Since if V> € Af( ( D + BB m )\s n ), 
then 

£>V> = 0; B 9 xl> = 0 

Then it is easy to see that (156) has the following solution 

x(t) = (a cosa; n t + 6sina; n <)V’ 
which is obviously not strongly stable. 

3. In particular, if we let f(x) = x, then we obtain that the necessary and 
sufficient condition for the linear system 

x(t) + Dx{t) + BB*x(t) + Ax(t ) = 0 

to be strongly stable is Assumption 3 holds (or Assumption St holds if all w n ’s 
are distinct). . 

4. Assumption 1 not only plays the role of assuring the existence of solution of 
(156), but also contributes to the strong stability of (156). In other words, to 
guarantee the strong stability of (156), it is important to require /(•) being 
continuous at least in the neighborhood of x = 0. In the following example, 
which is although of single-DOF, we can see since /(•) is not continuous at 
x = 0, x(t) does not always go to 0 as t — * oo, even Assumptions 2, 3 are 
satisfied. 


EXAMPLE: Consider the governing equation of a spring-mass system with 
Coulomb damping 

x(t) + /isgn(x(f)) + ulx{i) = 0 

in which H = JR}. 
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By considering x(t) > 0 and i(t) < 0 separately, we can obtain the solution 

r 

u> 0 x(f) + ;£sgn(x(f)) = u> 0 x(0)cosu>o< + x(0)sinu> o t 

+H sgn(x(t))cosu;ot/wo 

< ( 157 ) 

x(t) = — u>ox(0) sin tvot + x(0) cos u>o t 

—fi sgn(x(t))sinu;ot/wo 

Therefore, in the (u> 0 x(t), x(<)) phase plane, the trajectory is governed by the 
following circle 

[u> 0 x(t) -|- — sgn(x(t))] 2 + [x(f)] 2 = [u; o x(0) + — sgn(x(t))] 2 + [x(0)] 2 
u>o 

When the representative point (u>ox(<),x(f)) is in the upper half plane, the tra- 
jectories consist of a series of circular arcs with center (—fji/u> o,0) and radius de- 
pending upon the initial data or the state when the representative point enters the 
upper half plane from the lower half plane. Similarly, when the representative point 
(u;ox(<),x(t)) is in the lower half plane, its trajectories consist of a series of circular 
arcs with center o, 0), see Figure 2. 

If the initial data satisfies |u; o x(0)| > h/uj 0 , then, starting from the initial point 
(a> o x(0),x(0)), the representative point moves clockwise along various circular arcs 
in the upper and lower plane. The process stops when the representative point 
intersects the u>ox(t) axis between —h/uq and fi/u>o. The motion ceases at such a 
point because the maximum possible friction force exceeds the force in the spring, 
i.e. 

H > u>ox(<) 

Similarly, if the initial data is such that 

x(0) = 0; |x(0)| < fi/ul 
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then there is no motion because the spring force u;qx( 0) cannot overcome the friction 
force /x. 

Therefore, from this example we can see, since 

f(x) = /x sgn(x) 

which is not continuous at x = 0, the system energy does not generally go to zero 
as t — ► oo. After the mass stops moving, its strain energy 1/2 w£x 2 (f) is generally 
positive. 

PROOF OF THEOREM 13: First of all, Assumptions 2 and 3 imply /( 0) = 0. 
In fact, let 

fl (^-li " ' ) ^m) 

/(*) = i 

^ fm{ x ' ' ' >®m) 

Suppose some /,(0,---,0) ^ 0 for some j. Without loss of generality, suppose 
/i(0, • ■ • , 0) > 0. By the continuity of / at x = 0, there exist 6 > 0, e > 0 such that 

/i(xj,0, • • • ,0) > /i(0, • • • ,0) — e > 0, for |xj| < 8 

Then, 

M-6/2, 0, • • • , 0)(- 6/2) < -W, (0, • • • , 0) - e] < 0 

However, from Assumption 2 we already know 

fi ( — <5/2, 0, • • • , 0)(— 6/2) > 0 

This contradiction implies /i(0, • • • , 0) = 0. Similarly, we can show 

/i(0, ' - * ,0) = 0, i = 2, • • ■ , m 
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Next, for any (x(0),x(0)) € V(A) <g> V{D), by virtue of Assumption 1, the 
existence and uniqueness of solution of (156) is immediate. 

Since 

^ = -[Di(t),i(f)l - ( HB-x(t)),B-i(t)} < 0 
E(t ) is monotone decreasing and is lower bounded (by 0). Therefore, there exists 
E( oo) > 0 such that 

lim E(t) = E{ oo) 

t— »oo v 

It suffices for us to show that E( oo) = 0. For this purpose, let us assume that 
E( oo) > 0. Then there exists energy-preserving steady state motion, denoted by 
(x,(t),x,(<)), with corresponding energy 

E.(t) = 1/2 ([Ax,(<),x.(0] + ||x.(t)|| 2 ) = E(oo) 

Then, 

at 

= o 


which immediately implies 

( Dx,(t) = 0 
[ B m x,(t) = 0 

by Assumption 2 and /( 0) = 0. However, (158) indicates that x,(f) solves 


(158) 


x,(t) -I- Ax,(t) = 0 


or 

OO 

x,(f) = ^2(a n cosu n t + b n sinu n t)rj} n 

n= 1 
oo 

x,(<) = 5^u> n (6 n cosu; n < - a„sinw n <)V> n 

n=l 
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Since we do not intend to lose the generality by excluding the repeated u>„ case, 
we use {V’fc,} to denote the eigenvectors corresponding to the repeated natural fre- 
quency u> k . Then, x t (t) can be rewritten as 

OO 

X,(t) = cos u> k t(^T h, 'J>k , ) 

k—1 ki 

OO 

- sin u k t(%2 a ki rj> k , ) 

k=i ^ 


Then, (158) is equivalent to 

OO 

J2 u k cos u k t(D + BB*)(52b ki ip ki ) 

fc=l ki 

-^2u> k sin u k t(D + BB*)(%2 a*, V’t. ) = 0 (159) 

k=. 1 ki 

Therefore, it must be the case that 


(D + BB-)(Z„,b k M = 0 

. (C + BB-KE*, <■*,*.) = 0; t = l,2, -- 

Then, by Assumption 5, we know that 

^ b k ,tpk, = o 

i Sfcj a k,^ki = 0; k = 1,2, •• • 
Therefore, the orthogonality of {V’i,} implies 


(160) 


(161) 


a kx = 0; b kt = 0, for all fc,-, and k = 1,2,--* 


i.e., x,(t) — 0 or E( oo) = E,(t) = 0. 

From the proof we can easily see that if Assumption S is replaced by 

Assumption 3 a : [DrJ> n ,iJ> n ] >0, n = 1,2, • • • , 


□ 
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then Assumption 2 can be replaced by a weaker version, i.e. 

Assumption 2 a : [/(x),x] Rm > 0 for x € JFT 1 . 

Under Assumption 8a, the uncontrolled system 

x(t) + Dx(t ) + Ax(t ) = 0 

is itself strongly stable. By using active damping u(t) = — f(B*x(t)), we are en- 
hancing the system stability. 


7.4 Mode Excitation by Active Damping 

In this section, we study the following problem: For an undamped linear oscillation 
system, 

[ x(<) + Ax(i) = 0 


( x(0) = x 0 x(0) = io 

If x(0), x(0) are linear combinations of finite number of modes, say, 

K 

x 0 = ^2a nk rpn k 
*= 1 
K 

io = 

k=l 

then the corresponding solution x(t) always stays on these modes and is of the form 

K 

*(t) = E a » *(t)&n* 

k=l 
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where a njj (t) satisfies 

a n*(0+ w n*°n*(0 = 0 
' a n*(0) = OCn k 

. <(0) = ft.. 

However, with active damping, the system response no longer stays on the initial 
modes. For simplicity of notations, suppose the damped system starts from the first 
mode, i.e., 

x(0) = a!(0)V>i, x(0) = 0 

then the active damping will excite other modes. How many more modes are excited? 
In what magnitudes? What is the behavior of those excited higher order modes? 
These are some of the questions we will answer in the following analysis. The idea 
is to solve the feedback control effort z(t) = a finite dimensional time 

function, without solving the whole system - a PDE. 

First, we write the damped system response in its mode decomposition form 

*(0 = J2 a « Wn 

n=l 

where the mode responses a n (t), n = 1, 2, • • • solve 

' a n (t) + u;£a n (<) = —4> n (L)/b 2 z(t) 

* MOHrnCO), d|(0) = 0 (162) 

a n (0) = 0, a„(0) = 0, n > 2 

where z(t) = f(u(t,L)), the rate feedback control. 
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Recall that (151) can be written as 


w(t) = Aw(t) + Bz(t) 


where 



The semigroup T(<) generated by A is given by 


( C{t ) S(t) \ 


where C(t) : H — ► H, S(t ) : H — ► V{A l12 ) are cosine and sine operators defined 

by 


C(t)l = Y COSU> n t[l,V’n]V’n, 


n=l 

oo 


_ . . « > sin r / i i 

S(t)y = [y,^n]V»n, 


n=J U * 


x e H 
yeH 


Therefore the response (x(t), x(t)) is given by 


I(,) = T(i) I(0) 

i(f) / V ± 1 °) 


j + jf T(< — t)Bz(t ) di 

from which we obtain, noticing that i(0) = 0 and x(0) = <zi(0)^i, 

x(t) = C(*)x(0) + / S(t - t)(-Bz{t))<It 

J 0 

f 1 ^ sinu>„(* - r) <f> n (L ) 


ft w 

= ai(0)cosu>it tpi — I 5Z 

Jo n=l 


U>n 


fc 2 


^ n z(r)dT 


Then, separating the boundary component in x(t) gives 


1 /* ^ sinu> n (< — r) ,, . 

u(t,L) = a a (0)^i(L)cosu>it - r? / Y <t>l(L)z(T)d T 

° J ° n=l UJn 


(163) 
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Differentiating both sides with respect to t gives 


ti(f, L) = — ai(0)<^i(LVi sinwjt — f K\(t — t)z(t)<1t 

Jo 


(164) 


where 

Ki{t) = ^ Yi 4>l(L) cos u n t 

From (154) it is obvious that the infinite series is absolutely convergent and the 
convergence is uniform in t. 

Therefore, from (164) we conclude that the feedback control z(t) is uniquely 
determined by the following nonlinear integral equation: 

z(t) = /[— ai(0)<£i(Z/)u>i sina>i< — f K\(t — t)z{t)<1t) (165) 

Jo 

If we can solve z{t) from (165), then we can find each mode response a n (t) from 
(162). We first study 

(1) Linear Damping Case (/(x) = kx) We introduce the notation 


p/ \ 1 ^n(L)^>n(0 

G ( s - r ) = iiL 


b 2 s 2 + W 2 

n= 1 1 n 


0 < r < L, s € C 


And, later on, we will often use the capital letters to denote the Laplace transforms 
of the functions denoted by the corresponding lower case letters, such as 


Z(s) = £[*(01, A.C) = £M01 


etc. 


By performing Laplace transforms on both sides of (165), one cam obtain 

—kai(0)u*<f>i(L) 


Z(s) = 


(s 2 + u>?)(l + ksG(s,L)) 
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Furthermore, from (162), 


M s ) = 

A n {s) = 


one can obtain the mode responses 


*(1 + *•/**££> $$) + ktf(L)/l> 2 

° l(0) (s*+u*)(l + ksG(s,L)) 

fcai(0)u;?<^i(L)^n(^)/ft 2 
(s 2 + W?)(s 2 + u> 2 )(l + ksG(s, L)) 


n 


2,3, 


(166) 


From the definition of G(s, L), we can realize that {±tu; n , n = 1, 2, • • •} are not 
the poles of A n (s), n = 1,2, ••• due to cancellation. In fact, all A n (s)’s have the 
same poles and they are simply the roots of the equation 

l + ksG{s,L) = 0 (167) 


(2) Saturating Nonlinear Damping Case. First of all, the saturation func- 


tion /(*) can be written as 


f(x) = kx — tp(x) 


where ift(x) is the following dead-zone function 




0 |*| < M/k 

i 

k(x — sgn(x)M/fc) |*| > M/k 


Then from (164), we have 


z(t) = f(u(t,L)) 

= kii(t, L) — L)) 

= k[—ai(Q)<f>i(L)u>i sinu>i< — f Ki(t - t)z(t)<1t]- 

JO 

Then, performing Laplace transform on both sides and letting 


tf(s) = £[r/>(u(f,L))] 
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give 


+ «*?)*(*) 

K) (s 2 -fu> 2 )(l + *sG(s,.L)) 

Similarly as in linear damping case, we obtain 

... j(i + i./* , K^ £&) + *«£)/* *(£)/»»*(,) 

l ) ° lW (»’+«?)(l + JbG(*,£)) (3 ! + u;)(1 + *:8G(s,I)) 

. , . t»,(D)u. 1 V 1 (£)^(£)/y , A.ffl/WW nBS , 

" { ‘ ~ (s2+wf)(» J +«;)(i + ifc»G(i,£)) (» ! +w;)(i + *»G(«.£)) v 1 

n = 2, 3, • • • 


Comparing (168) with (166), we can realize that each mode response .A n (s) is 
simply the linear damping mode response A n (s) plus a correction term C„(s), with 


C n (s) = ML)/b 2 


*(f) 

(s 2 + u> 2 )(l + ksG(s, L)) 


Furthermore, all the poles of each correction term are again the roots of (167), 
because we can show that 'P(s) is analytic on the complex plane. 

In fact, from Theorem IS , we can see that under saturation type nonlinear damp- 
ing, the damped system is strongly stable. Hence, 


l/6|tt(/,£)| < ||x(f)|| < \p.E{t)] l l 2 — * 0, as t — *• oo 


That is, for fixed linear range slope k , and saturation level M , there exists T = 
T(k, M) < oo such that 


\ii(t,L)\< M/k, for t >T 


Therefore, ^(s) can be written as 

*(») = CMUt.m = f T <Hi.(t,L))e-“dt 

Jo 
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Since T is finite and rj>(u(t, L )) is bounded and continuous with respect to t, we thus 
conclude that 'P(s) is analytic. 

From above analysis we can see that in both linear and nonlinear damping cases, 
all /4„(s), n = 1, 2, • • • have the same poles, which are simply the roots of (167). 
Therefore, we now introduce 


Definition 1 (167) is called the Characteristi£ Equation of the distributed param- 
eter system (150). 

To study the behavior of the mode responses a„(<), n = 1, 2, • • •, we need to 
study the characteristic equation in detail, in particular its root locus. 


1. The closed form expression of G(s,L ) is given by 

W s 2 cosh \/~ L "f cos \J ~~ L 

i/g( S ,l)= s 2 +-^(-) —^7—7 jr~ 


Hence, the characteristic equation (167) is reduced to 

t 2 y/s 2 -+ cosh + cos 
>/2 sinh — sin 


(169) 


(170) 


In particular, we have 

£ ^ = b'G(0,L) = L 3 / 3 

U>1 

n=l n 

The derivation of the closed form of G(s , L ) is given in Appendix A. 


2. G(s, L) is analytic on C\{±u*>„, n = 1,2, • • •}. {±tu; n , n = 1, 2, • • •} are the 
first order poles of G(s, L ). All the zeros of sG(s, L) are on the imaginary axis, 
except for one at infinity, i.e., they are {0, 00 , ±iz„, n = 2, 3, - • •} where 


2 „ — | — ► 0, as n — ► 00 
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In fact, the first statement is obvious. To show the second statement, it is 
easy to see that 0 and oo are two of the zeros of sG(s, L) from the definition 
of G(s,L). The other zeros, from the closed from expression of G(s, L), are 
simply the roots of 

sinh \j2sJaL — sin yfelfaL = 0 (171) 

However, the roots of the transcendental equation 

sinh z = sin z 


are given by {(1 ± t)xfc, k = 0, 1, 2, • • •} where {xj,, k = 0, 1, 2, • ■ •} are the 
roots of 


tan x = tanh x 


Obviously, we have 


x k — (k + l/4)ir| — ► 0, as k — ► 00 


Hence, from the relation z = \^2s / a.L, we know that the zeros of (171) are 
given by 

<ii ( ‘ 


or 


{±ia(x k /L) 2 , fc = 0,1,2, ••• } 

Therefore, letting z n = a(x n _ 2 /I') 2 , n = 2, 3, • • •, we realize that the zeros of 
sG(s,L) are {±i2 nj n = 2,3, •••}, in addition to {0, 00 }. Furthermore, from 
the definition of G(s, L), it is not difficult to see that 


u> n _i < z n <! u? n , n — 2, 3, • • • 
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3. For 0 < it < oo, the characteristic equation (167) has a sequence of roots 
{s„,3^, n = 1,2, • • ■} such that 

3t[-s n ] <0; n = 1,2, ••• 

and 

R[in) * 0, |^[s»i] ^nl * 6, Ti ► OO 

In fact, first of all, the characteristic equation (167) cannot possibly have any 
roots on the imaginary axis, because ksG(s , L ) will be purely imaginary on 
the imaginary axis excluding {±iu; n , n = 1,2,- • •}. 

To see the rest of the claim, it is sufficient to plot the root locus of (167). 

As usual, the root loci start from the poles of sG(s,L), i.e. n = 

1,2, •••} and end at the zeros of sG(s,L), i.e., {0, oo, ±tz„,n = 2, 3, •••}. 
The angles of departures from the poles and the angles of arrivals at the zeros 
are all 180 deg with respect to the positive real axis. Therefore, we can realize 
that the root loci stay in the left half plane, since it cannot cross the imaginary 
axis. 

Since one of the zeros of sG(s, L) is at infinity, we know that one branch of the 
root loci extends to infinity along an asymptote. The angle of the asymptote 
to the positive real axis can be computed to be 180 deg. Furtheremore, in 
general, the root locus on the real axis always lies in a section of the real axis 
to the left of an odd number of poles and zeros (in our case, there is exactly 
one zero on the real axis, which is s = 0). Therefore, the root locus occupies 
the whole negative real axis. 
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Of course, the root locus is symmetric with respect to the real axis, since 


1 + ksG(s,L ) = 1 + ksG(s,L) 


The root locus of (167) is shown in Figure 3. 

To verify the statements on s„, the roots of the characteristic equation, we let 
w = \J2 IJaL. Then (167) can be written as 

,1 2 kL. sinhw — sinu? .6., 1 

( — I ) 1- (-) — = 0 

L aw 2 2 + cosh w + cos w aw 

Since the roots w n — yj2s n /aL are unbounded as n — ► oo, we can conclude 
that 

|u > n+ 2 - (1 + i)(n + 1/4 )t| — ► 0, as n — ► oo 

Hence 


s n = a / 2(^) 2 

L Pn 


= «*>„( 


«>n/(l +0 x2 
fin ’ 


Therefore, we can show that 


£[s„] — ► 0, |9[,J - w n | — ► 0, as n * oo 


Based upon the above study of the characteristic equation and its root locus, we 
can obtain the following conclusions: 

(/) In linear damping case, we have found the Laplace transform of each mode 
response A„(s), n = 1,2, • • •, supposing the beam starts from the first mode. Active 
damping excites all other modes. The active damping process is as follows: the 


118 



- ksG(s , L) = 0, for 



system energy initially possessed by only the first mode, 1/2 u>j<zJ( 0), is first partially 
shifted to all other modes, and then the energy acquired by each mode is gradually 
absorbed by the active damper. Quantitatively, the excited mode response a n (t), 
n > 2, is equivalent to the output of a stable system with the following transfer 
function 


H n (s) = 


u$ML)/P 


l + ksG(s,L)) 

corresponding to an initial impulse excitation with magnitude fcai(0)<£ n (.L). Each 
of the excited mode responses is proportional to the feedback gain k. 

The excited mode response in time domain can be found through inverse Laplace 
transform. Let T be a closed contour consisting of the imaginary axis and the 
semicircle with infinite radius on the left half plane. Then for n > 2, 


where 


a n(Q 



e‘*ds 


^J r A„(s)*«ds 

J^[Res, =ik {>l n (s)e* < } + ReSaai^iM^e*}] 

k= 1 


Y, p{ n) exp(3fc[s*]t) cos(9 , [s*]< + 0[ n) ) 
k = i 


( 172 ) 


= |Res..«{A.(»)}| 

^[ n) = Arg[Res, =tl { A n (s))] 


From (172) we cam see that each mode response a n (t) involves all frequencies, 
which are slightly smaller than the corresponding natural frequencies. As frequency 
goes higher, the damping effect of the active damper becomes weaker, because 
3?[sfc] — ► 0 as k — * oo. 
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(ff) In saturation type nonlinear damping case, we have shown that A„(s) is 
equivalent to the sum of the A n (s) corresponding to linear damping case and a 
correction term. Furthermore, all the correction terms have the same poles as in 
linear damping case, i.e., the poles are the roots of the characteristic equation. In 
other words, the >4„(s) in both linear and nonlinear damping cases have the same 
poles, but with different residues. 

In saturation type nonlinear damping case, the mode response can be similarly 
obtained (attention needs to be paid on the restriction t T(k,M) to guarantee 
the integrand vanishing on the infinite semicircle). 

a„(t) = £ p[ n) exp(R[s*]<) cos(3[s;t]f + 0[ n) ) 

with similarly defined as in linear damping case. 

Therefore, the existence of saturation in active damping only causes variation of 
amplitudes and phases, and does not result in any qualitative change in terms of 
mode responses. In particular, the component frequencies remain the same and no 
extra frequency is introduced. These confirm that saturation type nonlinearity is 
amplitude sensitive but frequency insensitive. 

(M) From the root locus ( Figure S ), we can see that when the feedback gain k is 
chosen such that Si and ij are located to the left of PQ , the base frequency oscilla- 
tion decays much faster than higher frequency oscillations. Then higher frequency 
oscillations soon become dominant in this case. 

When k is sufficiently large such that the first pair of roots are on the real 
axis but to the left of PQ, base frequency oscillation is then completely eliminated. 
However, it is impossible to eliminate any of the higher order frequency oscillations 
with only one actuator. In other words, active damping effect is significant only to 
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the base frequency oscillation. 

(IV) From (172) we can see that for any n > 2, there does not exist a > 0 such 
that 

|a n (f)| ^ const. e~° x 

That is to say, the total energy E(t) of the system does not possess exponential 
decay. 

However, if we use modal truncation as an approximation, the finite dimensional 
truncated system energy does decay exponentially in both linear and saturation non- 
linear damping cases. This is a difference between the original infinite dimensional 
model and the finite dimensional truncated model. 

Let us condiser the saturation type nonlinear damping case. Let the finite mode 
approximation of x(t ) be 

N 

*(0 ~ Yi “"(OV’n 

n=l 

Then {a n (t), n = 1, • ■ • , N] satisfy 

a„(<) + <t> n (L)ltfz N (t) + u 2 n a n (t) = 0, n = 1, 2, • • • , N (173) 

where 

ZN(t) = f(Y2 a n (t)<f>n(L)) 

n= 1 

By introducing the following notations 

fi = Diagonal {u>j , , • • • , } 

f oi(t) 

XN(t) = : 

L O^(f) 
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Bn 


' ML) ' 

4>n{L) ; 


(173) can be recast into the form 

x N (t) + i B N f(B- N x N (t )) + Slx N (t) = 0 


*7v(0) = 

< «,( 0) ' 

, Xjv(0) = 

' 4.(0) ' 


^ QAt(O) y 


^ ajv(0) y 


Applying Theorem IS , we can conclude that for any initial data 


||*n(<)II — * 0, as < — » oo 

Then, through similar steps as in obtaining (164) we can obtain 

1 /* 

a n (t) = — a„(0)u; n sin u n t + a n (0) cosu>„t - j cosa; n (t - r)<f> n {L)z N {T)dT 

n = 1,2, — , iV 


and, further that 

2A/(0 = fCT. hn{t)<t>n(L)) 

n= 1 

N 

= fc[52(®»(0) cos Unt ~ a n(0)u>„ sinu> n t)^„(£) 

n=l 

—— f 52 <t>l(L)cosu n (t - T)z N (r)dT ] - d n (t)4> n (L)) 

*> 2 Jo £1 „=i 


Therefore, 


Zn(s) = 


1 


1 + ksGs(s, L) 


£ ^(0). a.(oy (I) _ 
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MO 


where 


-i— [o n (0)» + a„(0) - ML)/b 2 Z N (s)\ 

s * uj* 

Pn (») ML) 

+ u>?)(l + ksGsi^L)) b 2 (s 2 + u>l)(l + ksGN{s,L)) 


Gn{s,L ) 


b 2 


f ML) 

nt^+^ 


<M0 = CM'ZMVML))} 

n= 1 


and each P n {s) is a polynomial of s of order 2 TV — 1, for n = 1, 2, • • • , TV. 
Let 


<7 = mini <n<jv { — &[s„] } > 0 


where {s„, s^, n = 1, 2, • • • , TV} are the roots of the following rational characteristic 
equation 

1 ■+■ IcsGn(s, L) = 0 

Then, through inverse Laplace transform we can see that 


MO I < M n e-"‘ 
MO I < M ^~ ot 


for some M n > 0. And hence, 

E„(t) = 1/2 £>M(< ) + <£(()) < AToe- 2 ", t > T(k, M) 

n=l 

for some Kq > 0. Therefore, we can find K > 0 such that 

E N {t) < Ke - 2ct , t > 0 

i.e., the energy of the modal truncation model decays exponentially. 
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7.5 Sensor Noise Effect 


Since in most cases, measurements axe corrupted by sensor noise, we consider the 
effect of sensor noise in this section. We axe particularly interested in the sensor 
noise effect on the steady state motion of the offset antenna. For simplicity, we 
assume that the sensor noise is a white noise process n(t) with constant spectral 
density 1 over the whole frequency range. Then the sensor output is given by 

u(0 = B m x(t) — <rn(t) 

= u(t,L) — an(t) 

If we use direct connection, the feedback control is 

*(0 = f(ii(t,L)-<rn(t)) 

Again, we use the continuum model (150), which, in this case with sensor noise, 
can be rewritten as 


u(t, r) + a 2 u m '(t, r) = 0 

u(t, L) - fcV"(t, L) + /(«(*, L) - an(t)) = 0 

u(t, 0) = 0 (174) 

u'(<,0) = 0 

u"{t,L) = 0 

where f(x) can be either a linear function, f(x) = kx, or, a saturation function 
represented by 

f(x) = Mtan~ 1 (kx) 

for our convenience. 


125 



Linear Damping Case 

In linear damping case, things are rather trivial because we have the luxury of 
transfer function - we can find the transfer function between the beam tip position 
u(t,L) and the sensor noise input N(t). 

In fact, by setting the initial conditions to be zero, we obtain as before 


u (i, L ) = i / £ SmuJn ^ — l 'l<j> 2 n (L)z{T)dT 

b Jo n=l w n 

where, in this case, z(t) = k[u(t,L) — N(t)]. Then, through Laplace transform one 
can find 

U(s,L) kG(s,L) 

N(s) ~l + ksG(s,L) 

In particular, if N(t ) = crn(i), a white noise process, then we have the spectral 
density of u(t, L) in stationarity: 

a 2 k 2 \G(iuj,L )\ 2 


— 


|1 + kiwG(iu>, L)\ 2 

a 2 


— OO < U) < oo 


Then from the spectral density, we can further calaulate the stationary variance 
of u(t,L) by integration. Figure 4 is a typical plot of the spectral density 
from which, we can see that if the sensor noise involves only those frequencies, 
{ z n , n — 1,2, • • • }, i.e. the zeros of G{iz,L), such that 


N(t) = ^2 Cn cos z n t + d n sin z n t 

n 

then the beam tip is motionless at steady state, while the beam itself is vibrating 
under the sensor noise excitation. 

Instead of going further, we switch our attention to saturation type nonlinear 
damping case. 
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Figure 4- The spectral density of the stationary beam tip motion. 


Saturation Damping Case 

In order to cast (174) into Ito form, we first employ the following approximation 


f(u(t,L)-on(t) « X)) - L))on(t) + cr 2 /2/"(u(t, X)) 

- M (^) 2 [i + (^(U)) 2 ] 2 

so that the antenna motion equation in (174) is approximated by 


u (t,L) + f(u(t, L)) + <r J /2 /"(*(<. i)) - 6V'(I, £) 

= f\u(t,L))an(t) (175) 


Next, in order to simplify the diffusion coefficient, we divide both sides of (175) 
by l/(Mk)f'(u(t,L )), which is always positive, to obtain 


u{t , X) + MD(ku(t , X)) + [Jfcu(<, X)] 2 [u(f, X) - 6 2 u'"(f, X)] 

-tV"(t, X) = Mkon{t) (176) 


where 

D(x) = (1 + x 2 )tan _1 x — (o X) 2 ^ a: € iR 1 

Before we proceed, we pause on the curve of D{z) versus z. First, D(z) is an 
odd function, hence it is sufficient to study D{z) curve for positive z. We can easily 
have, when ok < 1, 

D'(z) > [1 — (cfc) 2 ] + 2z tan -1 z > 0 z > 0 

Therefore, D(z) is positive and monotone increasing on (0, oo) when ok < 1. When 
ok > 1, 

X>'(0) = 1 - {ok) 2 < 0 
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That is, there exists zq > 0 such that D(z) < 0 on (0, zo) and D(z) > 0 on (zo 5 o°)- 
The D(z) curves corresponding to a few ck values are plotted in Figure 5. 

Next, consider that 

[ku(t,L)) 2 [i 'i(t,L)-b 2 u"'(t,L)) 

= (ti.) 2 u((, L) - bV(t, £)] 

= (i6) , u(i,I)|i(l),*(i) + M')] 

= (kbf^-B'i(t) 

Hence, (176) can be rewritten as 

u(t,L ) + MD(kB"x(t)) + (kb) 2 E(t)B*x(t) 

—b 2 u'"(t,L) = Mkan{t) 

and, from which, (174) is approximated by 

£(<) + MBD(kB*i(t)) + ( kb)*E(t)BB m x(t ) 

+Ax(t) = MkaBn{t ) (177) 

In what follows, we are only interested in the steady state behavior of (177). 
After the system reaching stationaxity, the total energy E(t) fluctuates around a 
constant energy level and E(t) becomes a zero mean stationary stochastic process. 
Therefore, we can intuitively see that, after reaching stationarity, among the two 
damping terms, the nonlinear damping term MBD(kB“i(t)) becomes dominant 
over the other damping term, which has a zero mean random damping coefficient. 
Therefore, we neglect the damping term (kb) 2 E(t)B B*x(t) in (177) so that the 
stationary structure response is approximated by the following Ito type stochaLstic 
distributed parameter system: 

x(<) + MBD(kB m x(t)) + Ax(t) = MkaBn(t) (178) 
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Figure 5. The curves of D(ki) for various ok values. 



From (178), through similar procedures as in Section 4, we can obtain 

x(t) = f S(t — t)B[—MD(Icu(t,L)) + AfAr<rn(r)]dr 
'Jo 

where we have set the initial conditions to be zero since we are only concerned with 
the stationary behavior. 

By separating the boundary (scalar) component in x(t), we obtain from the 
above equation 

u(f , L) = M ( K 2 (t — r)[— D(fcu(r, Z,)) + k<rn{T)]dT 
Jo 

where 

K 2 {t) = 1/fc 2 ]T ■ sin u n t 

n=l Un 

From (154), we have 

Therefore, K 2 {t) is a fast converging series and we can reasonably use first N terms 
as an approximation of the infinite sum to obtain 

u(t,L) = ^ j V sin w n (f - T)[-D(ku(T,L)) + kcrn(T)]dT 

& J ° n=l W n 

By defining 

y n (t) = / 7? ^" smu; n (f - t)[— D( fcu(r, L)) + fc<rn(r)]dr 

Jo 0 U) n 

n = 2, 3, * * * 

we can realize that the beam tip (antenna) motion (u(t, L),u(t, L)) satisfies the 
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following system of stochastic differential equations: 


«(*,£) + gE?=,«£)]-D(M‘.£)) + of «(<,£) 


+SziM-“’W | ) = p!En=i OTMO 


Sn(0 + + oiJy„(() = g^J(i)t<rn(() 


(179) 


n = 2,3 , •••,A r 

We call (179) the ,/V-th order approximation of the beam tip motion. Notice that 
in (179), the damping coefficient 


N 


v E 


n=l 


because 


Y 2 EMr)ML) = \ 


n=l 


1 as N — * oo 

1 r = L 
0 0 <r<L 


In fact, it can be easily seen from 


= El 

n=l 


1 




= pE ML) ( Ur) 

6 ~1 V ML) 

In particular, the first order (N = 1) approximation of the beam tip motion is 
given by 

u(f, L) + \D{ku(t , L )) + u) 2 u(t, L) = A kan{t) (180) 

where A = M/b 2 <j>l(L). 
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When ka > 1, the deterministic version of (180) has a stable limit cycle ( Figure 

5 ). 

Stationary Statistics for Single-DOF Nonlinear Vibration 
In the rest of this section, we will find the approximate stationary probability 
density and varaince of (180). This is a central topic in nonlinear random vibration 
theory. Various methods have been reported by many authors. For surveys, see [25] 
[33]. 

The method we are going to use is called Method of Energy Approximation, 
developed by the author [72]. 

What we propose is to approximate a general nonlinear damping model 

x(t) + 7D(x,i) + u;ox(t) = aw(t) (181) 


by the following energy type nonlinear damping model 

x(<) + ( 182 > 


where p{E) is chosen in such a way that it minimizes 



[D(V2E/uo sin Vs V2 E 


COS 


V>) — p(E)y/2E cos rj>] 2 dip 


In the above minimization, E is considered to be fixed, because after reaching sta- 
tionary, the energy fluctuates around a constant level. In other words, at station- 
ary, the energy absorbed by the damping mechanism is statistically equivalent to 
the energy input due to the external noise. 

It turns out that p(E) is given by [72] 


p(E) = — f 1 D(y/^ /u 0 sin xj>,\/2E cosrp) cos ijtdip (183) 

ttV2E Jo 

The following are two facts supporting the above approximation. The first fact 
is [72] 
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Theorem 14 If D{x,y) is even with respect to x and odd with respect to y, then 
both the original model (181) and the corresponding modified model (182) have the 
same Krylov-Bogoliubov approximation, given by 


( da(t)/dt = —'y/LO 0 L(a(t)) 
\ d<p(t)/dt = 0 


where 


1 f 2ir 

L(a) = — / D (a sin xp,acoo cos ip) cos ip dip 

27t Jo 


The second fact is concerned with limit cycle. 

For the general nonlinear damping model (181) without noise input, i.e., 


x(t) + 7 D(x, x) + u>lx(t) = 0 


by Krylov-Bogoliubov approximation, we know that a limit cycle exists if and only 
if 

L{a) = 0 

has a positive solution a p . And in this case, the limit cycle has approximate ampli- 
tude a p and frequency wo- 

To find the stability of the limit cycle, let us define the deviation A a(t) from the 
limit cycle amplitude, i.e., 

a(f) = a p + Aa(t) 


Substituting into a(f) = —~f/u>oL(a(t)), we obtain 


dAa(t) 

dt 


—^/uoL(a p + Aa(t)) 
-lM^\*=a P )Aa{t) 
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Then, it becomes obvious that the limit cycle is stable if and only if 

dL(a) 


dt 


|o=o r ^ 0 


In summary, for a general nonlinear damping model, the existence, amplitude 
and stability of a limit cycle are determined by 

j L(a) = 0 
l > Oor <0 

Since both the original model and the corresponding modified model have the same 
L(a) function, we conclude that both of the nonlinear damping models have the 
same limit cycles and stability propertity. Of course, this analysis is based upon the 
Krylov-Bogoliubov approximation. 

For the modified model (182), its stationary probability density function is given 


by [72] 


2 f- 

P»(x,y) = Cexp[-— J 




fi(z)dz] 

1 /C = — /”exp[-4 /V(*)*M> 

UJq Jo <7* Jo 

In order to compute the function n(E) for our problem (180), we will need the 


following three identities 

fir/2 


fT { i 7r , . . 

/ tan _1 (fccosx) cosxdx = — r(\/l + Ar 2 — 1) (184) 

Jo 2k 

j * tan ~'[k cos x ) cos 3 xdx = + fc 2 — 1/2 + — - — ) (185) 

r ^t**.- * - - i /' /r + p ) ( i86 > 

Jo 1 + k 2 cos 2 x 2k* 


1 cos 2 x 


Using (184) - (186), we obtain 


n(E) = ~^== J XD(kV2E cos ij>) cos rpdxf) 


(ak) 7 


ME 


2 kE 2kEVT+2k*E 
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Next, we need to evaluate Jq fi(z)dz. For this purpose, we first calaulate 


Ij- 1ST + 2kzVl+ 2 Vz ]iZ = cr 2 ^ In(l + VI + 2k 2 E) + <7 2 fcln2 

= |(‘I + 2^)vTT2PS 

~ln(l + vT+2Pl) 

—k/2E + -^j-(2 In 2 — 8/3) 

3A: 

Then, we obtain 

"(A hyio^ )iZ = -^(4 + 2^)(8/9v^+2F1-1) 

+ -^(1 + -^j 2 )ln(l + y/\ + 2 k*E) + constant 

Therefore, the stationary probability density is given by 

p(z,y) = C[1 + 0 + fc 2 ( w?x 2 + y 2 )]° 

x exp{-/?[4 + k 2 (u\x 2 + y 2 )][8/9^/l + fc 2 (u;?x 2 + y 2 ) - 1]X187) 


where 


° = js + !^ 


j9 = 


1 


Xk(<rk ) 2 

and C is the normalizing constant, which is given by 

2x r°° 


1/C = X (1 + a;) 0 exp[-8/9/?(x 2 + 3)(x — 9/8)]da 


( 188 ) 


If we define 


/ n (a, $) = f x n (l + x)“ exp[-8/90(x 2 + 3)(x - 9/8)]ch 
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for n = 1,2, • • *, then 


c = ^/ r ' (*,/>) 


(189) 


It is not difficult to see that the stationary variance is given by 
P 2 _ 1 r hjo^) _ ,1 

E 2u\V [ h{*J) 1 

1 / 1 (a + 2,/?)-2/ 1 (tt+l,j9) 

2u\k 2 h(ot,P) 

Next, we study how the stationary density p(x,y ) is affected by the value of the 
parameter ak. 

For this purpose, we define 


((E) = y/l + 21 WE 

q( () = C(l+0°exp[-8/9^ 2 + 3)U-9/8)] (190) 


Obviously, 

p(x,y) = ?(\A + k^ufx 2 + y 2 ) ) (191) 

Since 

OOf+l 

sum = «'< i) = - (<7*n 

we realize that, when ak < 1, the stationary density p(x,y) attains maximum at 
the origin. For a typical plot of p(x, y), see Figure 6. When erk > 1, the maximum 
is not achieved at the origin, but rather, along a parabola with center located at the 
origin. A typical plot of p(x,y ) in this case is shown in Figure 7. Figure 8 is the 
density plot of p(x,y) corresponding to Figure 7, clearly indicating the altitude of 
the stationary density. 

We can actually find the equation of the parabola on which p(x,y) achieves 
maximum in the case ak> 1. Consider the equation 


,'({(£)) = o 
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i.e. <rk < 1 case. 



i.e. AAr > 1 case. 


1.5 



Figure 8. The density plot of p(x,y ) with er = 1, k = 3, and A = 1. 


which is equivalent to 


{ 3 + l/4(f 2 + £) = 3/2(u*) : 


(192) 


The solution of (192) is given by 


11 


£o(ot) = M3/2(cr*) J ) - IM3/2 (orfc) 1 )]- 1 - 1/12 


where 


= < 1^8 + 1/2 + ^ + < 1^8 + 1/2)211/2)1,3 


Therefore, the equation of the parabola is 


v'l + fcVi^ + y 2 ) = £o{<rk) 


or 

u?, 2 x 2 + y 2 = “ !] 

And, it is easy to verify that 

£o(ak) > 1 <=$■ ak > 1 

In Figure 9 , the variation of the stationary statistics Ex 2 as a function of k and 
A is plotted. It is obvious that the stationary variance increases more significantly 
in k direction than in A direction, for fixed observation noise intensity <r > 0. 

Among those parameters, a ■ the sensor noise intensity, and A - the maximum 
output of the actuator, are dependent upon the sensor and actuator used. The 
only parameter we can adjust is k in the linear range slope A k. Through the above 
analysis, we know that when k is large (ak > 1), the stationary variance (deviation 
from the equilibrium) is large. Therefore, to achieve small deviation from the equi- 
librium, a small k is desired. However, at transient, if k is small, the active damper 
no longer significantly enhance the stability. Therefore we have two options: 


141 




• preset a compromized value for k such that k « 1/<t; 


• let k vary in the following manner: 


k(t) = 


relatively large k\ at transient or when signal — noise ratio is large; 
relatively small &2 after transient or when signal — noise ratio is 
small. 


7.6 Conclusions 


A group of sufficient conditions for strong stabilizability is provided for general 
distributed parameter oscillation system, taking the actuator saturation into con- 
sideration. These are the weakest sufficient conditions obtained so far and it is 
found that the nature of internal damping is not crucial in guaranteeing the strong 
stabilizability. 

By extending the notions of Characteristic Equation and Root Locus to our dis- 
tributed parameter system, we studied the nature of active damping. Active damp- 
ing excites all other modes even if only one mode is involved initially. Saturation 
type nonlinear damping does not generate any qualitative change in terms of mode 
response. In particular, each component frequency remains the same and no extra 
frequency is introduced. Indeed, the Laplace transform of each mode response in 
both linear and nonlinear damping cases have the same poles, which are the roots 
of the characteristic equation. Through the root locus, it is also found that using 
one active damper can only significantly damp or eliminate the base frequency os- 
cillation and higher frequency oscillations become dominant in the response of the 
actively damped system. 

We studied the effect of sensor noise on the steady state response of the beam 
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tip. By calculating the approximate stationary probability density of the beam 
tip motion (u(f, L), u(t , L)), we have found that a large feedback gain (linear range 
slope) will result in a relatively large steady state deviation of the beam tip from the 
equilibrium due to sensor noise excitation. Therefore, in the design of the feedback 
gain, a compromise has to be made between satisfactory transient damping effect 
and reasonably small steady state deviation of the beam tip. 


144 




Chapter 8 


Conclusions and Some Open 
Questions 


8.1 Conclusions 

Frequency Response of Nonlinear Damping Model: Single DOF Case 


A method of computing the correlation function and the spectral density of 
nonlinear damping model is obtained. In the case that D(x,y) being of the form 
the spectral density is given by (40) in which ip v is given by (48). In 
the case of general D(x , y), the spectral density is given by (32) for which one needs 
to evaluate xj> v and m 2 ,o- V’v can be obtained from (31) or from (48) by first finding 
the corresponding /*(•) through the formula (45). Based upon (38), the approximate 
value of m 2 ,o is given by 


m2,o = 


1 m(l) 


u>o m(0) 

This paper also proposes a method to obtain the approximate explicit stationary 
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density for nonlinear damping model with general D(x,y). The idea is to replace 
D(x,y) in the exact model by the corresponding )y to obtain the modified 

model. The approximate explicit stationary density is given by (52) and (53). It is 
shown that both the exact model and the modified model have the same Krylov- 
Bogoliubov approximation and the same xp T , Vy 

In the Spacecraft Control Laboratory Experiment (SCOLE) program, the pri- 
mary control task is to rapidly slew or change the line-of-sight of an antenna attached 
to the space shuttle orbiter, and to settle or damp the structural vibrations to the 
degree required for precise pointing of the antenna. The objective will be to min- 
imize the time required to slew and settle, untill the antenna line-of-sight remains 
within the prescribed angle. From practical consideration, the maximum moment 
and force generating capability of the controllers on both the shuttle and the an- 
tenna beam/reflector are limited (maximum moment on both shuttle and antenna 
reflector is 10 4 ft-lb for each axis, maximum control on the reflector is 800 lb). 
Therefore, saturation type of control is inevitable. To avoid significant excitation 
of the beam while applying the slew control, the first harmonic of the slew control 
versus time should stay away from the resonant frequencies of the first a few modes 
of the antenna beam. This consideration is important in the design of the slew 
control. This paper provides an analytical frame of finding the spectral density of 
each mode, which is basically the amplitude of the system response corresponding 
to the sinusoidal input with each frequency u>. Specifically, the spectral density tells 
the designers where the resonant frequency is. 

Theorem S gives an interesting result concerning the modelling and identification 
of nonlinear internal damping in flexible space structures. In spite of the fundamen- 
tal importance of the damping term, the nature of internal damping has been little 
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known. In [8], the following nonlinear damping model is proposed 


x + 2(u n x + 7 x 2m |x| a x 2 " +1 |x|0 + a; 2 x = 0 


(193) 


The free response of (193) fits NASA flight data with great accuracy. However, if 

o 2 i i3 

we replace the nonlinear damping term in (193) by — )x where 

4 W 2 v/2 E , 

u(E ) = — 7 = / D( sin 0, y/2E cos 9) cos 6d6 

’ tV2EJo V w„ 

- JO-MJ + Sy 


zj2m-f q 

n 


in which 


9 


m + n + (a + /?)/2 

2r(m + 2±l)r(n + l + ^) 

7T T(q + 2) 


we obtain the following energy type nonlinear damping model 


x + 2fu > n x + 7 2 l° +a ( u l x2 + + W n* = 0 ( 194 ) 

By Theorem 3, both (193) and (194) have the same Krylov-Bogoliubov approxima- 
tion, that is to say, these two nonlinear damping models can not be distinguished 
based upon their free responses (for details, see [72]). 

Frequency Response of Nonlinear Damping Model: Multi-DOF Case 


A formula for computing the spectral density matrix of nonlinear damping model 
with n-DOF is presented, (77). The error of the formula (77) is of the order 0(7 2 ). 

It is not surprising that the spectral density matrix depends on the first order 
statistics R X i(0) and i? rv (0) for which we need the (first order) stationary probability 
density. However, to find the stationary density in general is itself a difficult task. 
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In the case of energy type damping, i.e. 

, ,x T KDx + y T MDy 
D(x,y) = fi( 


)<rcr T Dy 


the explicit stationary probability density can be obtained and is given by (90) and 
(91). In this case, the stationary density function is a function of energy Ed with 


Ed = 1/2(x t KDx + y T MDy ) 


and R xy ( 0) = 0, i? rr (0) can be obtained explicitly and is proportional to K *, see 

(95). 

As we know, in single-DOF case, x and x are always uncorrelated in stationary 
state. However, this luxiury does not extend to multi-DOF model. It is pointed 
out in this paper that even for linear multi-DOF model, x and x are in general not 
uncorrelated. A necessary and sufficient condition for uncorrelatedness is given by 
(82). 

The conclusions on the infinite dimensional model 


1. The shape of_ the frequency-response curve 

Let us consider any fixed mode, say the first mode. For single frequency 
excitation with frequency close to u\, i.e., 


uj = ui + ea 


the frequency-response equation is given by 

2 / 2 \ ( B ^ 1 /( 2 o > 1 )) 2 

M j ~ a 2 + [fo* + 7o/2(u>M(<t 2 ))*] 2 

The unique positive solution, denoted by g(cr\u)i,B' , tf)i), when plotted, is still 
bell-shaped. And there is no multi-peak phenomenon because g(-\u>i,B*xl>i) is 
an even function of a and is monotone decreasing as a 2 increases. 
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The frequency-response curve is given by 


I </(<7|u>i, B m \j> i) when u> = u>i + ea 

c lS^l when = °(i) 

or, generally, 


/>»(“>) 


g(cr\uj n ,B m xl> n ) 

< 




when u = w n + ta 


when |u> — u> n | = 0(1) 


However, for nonlinear stiffness (linear damping) problem, the shape of the 
frequency-response curve is more complex and qualitatively different. For 
example, for the following Duffing oscillator 


x + 2 cfii + J^x + eax 3 = £(<) 


the corresponding frequency-response equation is given by 

9 


One c an easily see that p(o) is no longer an even function of o'. In fact, in 
certain range of frequency, p(cr) is even multi-valued. As long as a ^ 0, the 
frequency-response curve is a backbone curve. 


In this case, there are jump phenomenon and caotic behavior which will be dis- 
cussed later. Of course, if a = 0, the model becomes linear and the frequency- 
response curve is single- valued and takes the shape we are familiar with. 


2. Comparison with linear damping problem 


149 


It is well known that for a single DOF linear oscillator 


x + 2e£ou>ii + u> 2 x = i cos u>t 
the stationary amplitude is given by its spectral density 


pK *) = 


(u 2 — u\) 2 + 4(oC 2 WjW 2 


When u> = u>\ + e< 7, i.e., 


u> 2 — u % = 2eau>i + e 2 cr 2 


we have, for the linear model, 

pIM = 




( ° + C 2^) 2 + (frwi + 

While the corresponding nonlinear counterpart is 

(*•*/( 2^i )) 2 


tfrM = 


<7 2 + 1^1 + 7o/2(a> 2 p^r)’] 2 


As we can see, 

Pn(u>) < Pl(u) for u; = u>i + e<r 

This is not surprising because inclusion of nonlinear damping makes the to- 
tal damping greater. This in turn results in smaller steady state response 
amplitude. 

When w, the excitation frequency, is away from u>i, we have 

Pl(v) = 


= < #f^ +o(£3 > 
= e ]in^|j +0(£2) 
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i.e., these two frequency-response curves are very close when the excitation 
frequency is away from the natural frequency. 

From above we can see, comparing with linear model, the frequency-response 
curve of nonlinear damping model has no qualitative change and the only 
difference is the frequency- response curve of nonlinear damping model is below 
the curve corresponding to the linear model, especially in the neighborhood of 
the natural frequency of that mode. 

3. Jump phenomenon and chaotic behavior 

For nonlinear stiffness problem such as the Duffing oscillator , the multival- 
uedness of the frequency-response curve due to the nonlinearity of stiffness 
has a significance from the physical point of view because it leads to jump 
phenomenon. To explain this, we imagine that an experiment is performed in 
which the amplitude of the excitation is held fixed, the frequency of the exci- 
tation, i.e. o, is slowly varied up and down through the natural frequency uq. 
We observe the amplitude of the harmonic response. If a starts from the left 
side of the peak and increases, the amplitude will jump from the peak value 
to the lower value. Conversely, if a starts from the right side of the peak and 
decreases, the response amplitude will jump from the lower value to a higher 
value. This jump phenomenon is due to the presence of nonlinearity. 

Then what value does the steady state amplitude take if the excitation fre- 
quency starts and stays at a point within the multivalued region? The answer 
is, it depends on the initial condition. In other words, if more than one steady 
states exist, the initial condition determine which steady state is physically 
realized by the system. 
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This chaotic behavior is exclusively possessed by nonlinear stiffness model. For 
nonlinear damping (linear stiffness) problem, the frequency-response curve is 
always single- valued. This means that the steady state response of a nonlinear 
damping system is independent of the initial conditions. And for nonlinear 
damping problem, all the points in the frequency-response curve are physically 
realizable and there is no jump phenomenon or chaotic behavior. This is 
one of the fundamental differences between nonlinear damping problems and 
nonlinear stiffness problems. 

4. Internal resonance and energy exchange between modes 

What is internal resonance? For multi-DOF nonlinear systems, an important 
case occurs whenever two or more natural frequencies are commensurable or 
nearly commensurable. Examples of near-commensurability Me 

U >2 ~ 2u>i , u >2 ~ 3u>i , 0^3 ~ U >2 ± U >1 , 

Lt > 3 ~ 2u>2 ± U >\ , U>4 « W3 ± U>2 ± 

Depending on the order of the nonlinearity in the system, these commensurable 
relationships of frequencies can cause the corresponding modes to be strongly 
coupled, and an internal resonance is said to exist. When an internal resonance 
exists in a free system, energy imparted initially to one of the modes involved 
in the internal resonance will be continuously exchanged among all the modes 
involved in the internal resonance. 

For example, we consider the motion of a mass m attached to a spring that 
is swinging in a vertical plane. If we let x(t) denote the stretch in the spring 
beyond its equilibrium and 0(f) denote the angle between the spring and the 
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vertical line, then the governing equations of the motion are 

, 0 sin 0(f) + 2x(f)0(f) _ ft 

m + TTW) ” 

x(f ) + k/mx(t) — (/ + x(t))6 2 — g cos 0 = 0 

where Jfc is a spring constant, l is the natural length of the spring, and g is the 
acceleration of gravity. 

The two natural frequencies are 

w, = (s/0 1 ' 2 

u> 2 = (k/m) 1 ^ 2 

Suppose l and m are chosen such that u>2 « wj. If one starts the motion 
when 6 = do ^ 0, by pulling the mass m down, one finds that the mass 
oscillates up and down first, and that then a pendulum-type component of 
motion develops and grows at the expense of the spring- type motion. After 
a while, the pendulum-type motion starts to decrease and the spring-type 
motion starts to grow. Thus the energy is transferred continuously back and 
forth between the two modes of oscillation. 

Whether commensurable or nearly commensurable frequencies can cause in- 
ternal resonance depends on the degree of the nonlinearity and the geometry 
of the system. 

For energy type nonlinear damping system, internal resonance never occurs 
for any commensurable or nearly commensurable frequencies. This fact can 
be seen from its Krylov-Bogoliubov approximation. By Krylov-Bogoliubov 
approximation, we know the energy possessed by the nth mode is given by 

£n(0)(l+79£ 9 (0)<)-£ 
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where E n (0) is the initial energy of the nth mode. Therefore, the energy of 
each mode is continuously absorbed by the damping mechanism and there is 
no exchange of energy between any two modes. This is another fundamental 
difference between nonlinear damping problems and nonlinear stiffness prob- 
lems. 


5. Steady state response to multi- frequency excitation and coupling effect 
Suppose the external excitation contains M frequencies, say they are 

uj n + ecr n n = l,2, 


Then the steady state response is dominated by these M modes. And the 
frequency- response equation for the nth mode is given by 

(/n£*V>n/(2u>„)) 2 


pIM = 




where G(d 2 ) solves 


G(S 2 ) = 70 / 2 E 

n=l 

and we used the notation d 2 = (<7j, 


(/ n g^ n /2) 2 
+ + G(<t2))2 J 


Comparing with the single frequency case, the stationary amplitude corre- 
sponding to multi-frequency excitation becomes smaller due to the coupling 
effect of nonlinear damping. However, qualitatively, there is no change in the 
shape of the frequency-response curve. Therefore we see that the coupling 
effect due to nonlinear damping is weak. This can also be seen by examining 
the Krylov-Bogoliubov approximation, which tells us, to certain accuracy, if 
the initial data involve only finite number of modes, nonlinear damping itself 
will not involve modes other than those involved initially. The free response 
will stay on those modes initially involved. 
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In one word, by using nonlinear damping, one can only quantitatively change 
the decay rate so that test data can be fitted. Nonlinear damping model with linear 
stiffness does not produce any peculiar behavior such as internal resonance, jump 
phenomenon or chaotic behavior. If the experiments on flexible space structures 
indicated the existence of any of these peculiar behaviors, nonlinear stiffness model 
becomes necessary. 


8.2 Some Open Questions 

Inspite of the above investigations, some questions remain unsolved. They are listed 
below: 

1. What is the spectral density of the following nonlinear stiffness model with 
linear damping 

x + 2£x + g(x) = <rn(t) (195) 

where £ > 0 is a small parameter, and g(-) is an odd function of x € St 1 . 

Even more challenging is the same question without assuming £ being small. 
A handy example is the so-called Duffing oscillator 

x + Tfi 4- x 3 = an(t) ( q > 0) (196) 

R. N. Iyengar [45] studied this problem by first enhancing the dimension 
(DOF), then using Equivalent Linearization Mehtod(ELM) to obtain an ap- 
proximating linear system of equations, thus reducing the original single DOF 
nonlinear stiffness model to a two-DOF linear model, for which the spectral 
density become trivial. To illustrate the idea, let us consider the above Duffing 
oscilla'or (196). 
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Let z — x 3 , then 


i = 3x 2 x 

z — 6xx 2 + 3x 2 (<rn(t) — x 3 — tjx) 

= 6xx 2 + 3 ax 2 n(t) — 3 zx 2 — rji 

Therefore, we obtain the second equation 

z + r)z + 3 x 2 z — 6 x 2 x — 3a x 2 n(t) (197) 


Then by applying Equivalent Linearization Method to (196) and (197), one 
obtains the corresponding two-DOF linear model 



where P, Q, C\ are certain positive constants. 



(198) 


The spectral density obtained by this approach shows two peaks reflecting the 
existence of subharmonics in the system. The secondary resonance occurs at 
about three times the primary resonance frequency. Notice that the resulted 
linear model has non-symmetric stiffness matrix. 


Only Equivalent Linearization Methoditself cannot account for the existence of 
higher harmonics which is one of the very important nonlinear phenomenon. 
Applying ELM to a single DOF nonlinear oscillator leads to a single DOF 
linear oscillator which can oscillate at only one frequency. This viewpoint 
hints the desirability of increasing the DOF of the equivalent linear system, so 
as to allow more than one natural frequency to exist. 
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It is worthwhile to mention that the stationary probability density of nonlinear 
stiffness model (195) has been found, see, for example [25], which is given by 

P$(x,y) = Co exp[— 2£/ a 2 (G(x) + y 2 /2)j 

where Co is the normalizing constant and 

G(x ) = I g(u)du 
Jo 

2. In Section S of Chapter 8, to find the approximate explicit stationary probabil- 
ity density of general nonlinear damping model in single DOF case, what the 
author proposed is to replace the damping term D(x, y) by the corresponding 
energy type damping )y, where p(E) is given by 

p(E) = — 7= / 1 D{y/2E/ix>o sin xl>,V2Ecosi}>)cosrl>di/> (199) 

7 r\/2 E Jo 

And we have seen that both 

x -I- £D(x, x) + WqX = <rn(<) 

and 

X + (f*(— — 2 ) x + U 0 X = <rn(t) 

have the same Krylov-Bogoliubov approximation and the same rp x and tp v , 
which ju-e defined in (31). What we do not know now is the difference (in 
an appropriate sense) between the two stationary probability densities corre- 
sponding to the above nonlinear damping models. And what are the differences 
between the stationary (first order) statistics, such as i?*x(0) and R in { 0)? 

In this regard, I make the following conjecture: for the general nonlinear 
damping model with single DOF 

x + D(x, x) + u>qX = an(t ) 
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the stationary probability density can be written in the following form 

2u> 0 ff D(u,v)v 


p,(i,y) = C 0 exp[ -ff 

T(T*J J (u t v) 


-dudv] 


(200) 


(u,v)€S(x,y)U>oU 2 + V 2 

where S(x,y) is a (x, y)-dependent closed area. 

As we already know, the stationary probability density for energy type non- 
linear damping model is given by 

f E 

p,(x,y) = Co exp[— 2/ < 7 2 / fi(z)dz) 

Jo 

If we use (199), the relation between p(E) and D(x,x), then we obtain 

rE 

P*(x,y) = Coexp[-2/<7 2 / (i(z)dz ] 

Jo 

2 [E r2* y/2z r— 

= C 0 exp[ / / D( sint/>, v2zcosi/>) 

7T<7 2 JO JO O?o 


2 f* f*« ^y/Tz 
<jJq 

X COS t{>dtj>dy/2z] 

2 /nA'o^+V 2 f 2 * , , 2 


2 / , V u 'o x+w 2 

= C 0 exp / / 1/p" 

JO JO 


xZ)( 

2u?o 


/.\ _2 


O>0 


, p cos V>)p 2 cos tpdipdp] 


= Co exp[ \[f . 

ira 2 JJ u*u*+v3<uZx 


D(u, v)t> 

U$U 2 + tf 2 


dudu] 


In this particular case, S(x, y) is a (x,y)-dependent parabola with center at 
the origin and axes (w£x 2 + y 2 ) 1/,2 /w o and (w^x 2 + y 2 ) 1 / 2 . 

3. If the self-adjoint operator A is positive definite and has compact resolvent, 
with 0 € p(A), then the eigenfunctions of Ao 


A def 

Ao = 


0 1 

-A 0 
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given by 





K 

lU>n<i>n 

form an orthogonaJ basis in 

In fact, first it is easy to verify that {$+, tf'jn 
sequence. To prove its completeness, let 


; n = 1,2,---} 


1,2, • ■ •} is an orthogonal 


and 



eH®H 


w _L span{$+,$„; n = l,2,---} 
. Or equivalently, we have 

K *£]je; = <4[x, <f> n \ - tu; n [y, <f> n ] = 0 
[ty, = W^[x, <j > n I + tw n [y, <f>n I = o 
which immediately gives 


[y><£n] = 0; [x,^ n ] = 0, n = l,2, ••• 

By the assumptions upon A we know {<f> n \ n = 1,2, •••} is an orthonormal 
basis in therefore, x = 0, y = 0, i.e., w = 0. Then the completeness is 
justified for the non-damped case. 

Now the open question is, with damping term D , are the eigenfunctions of 

AH ( ° ' 

\ -A -D 

complete in H ® Til 
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